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Abstract. 

Lattice field theory methods, usually associated with non-perturbative studies of 
quantum chromodynamics, are becoming increasingly common in the calculation 
of ground-state and thermal properties of strongly interacting non-relativistic 
few- and many-body systems, blurring the interfaces between condensed matter, 
atomic and low-energy nuclear physics. While some of these techniques have been 
in use in the area of condensed matter physics for a long time, others, such as 
hybrid Monte Carlo and improved effective actions, have only recently found their 
way across areas. 

With this topical review, we aim to provide a modest overview and a status 
update on a few notable recent developments. For the sake of brevity we focus 
on zero-temperature, non-relativistic problems. After a short introduction, we 
lay out some general considerations and proceed to discuss sampling algorithms, 
observables, and systematic effects. We show selected results on ground- and 
excited-state properties of fermions in the limit of unitarity. Appendices include 
technical details on group theory and selected derivations. 
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1. Introduction 

The term 'lattice methods' is commonly used to refer to a set of computational 
physics tools whose unifying idea consists in describing the spacetime continuum as a 
(generally hypercubic) mesh, of finite volume V and finite spacing b. Quantum fields 
are then defined on such a spacetime grid and the many-body problem thus becomes 
a problem with a finite number of degrees of freedom. As a matter of convenience, the 
corresponding partition function is typically represented in terms of Euclidean path 
integrals. 

The fundamental reason for employing lattice methods to treat a given problem 
is that in many cases of interest the resulting path integral can be evaluated 
stochastically using computers, in a fashion that is both efficient and in principle 
exact. Indeed, assuming the problem is such that a proper probability measure 
can be defined, the stochastic calculation of the lattice path integral will suffer only 
from statistical uncertainties, which can be precisely quantified. Furthermore, errors 
associated with the use of finite volumes V and finite lattice spacings b, are also 
controlled and can be systematically removed by using larger V, smaller 6, and/or 
improved actions. Thus, lattice methods, when applicable, surpass perturbative as 
well as mean-field approaches in the quest for the reliable quantitative characterization 
of quantum many-body systems. This is true in particular if those systems are close 
to a phase transition, or if they are in any way strongly correlated, by which we mean 
that the low-energy excitations (and thus the long-distance properties) of the system 
are dramatically different from those of a system with the same degrees of freedom 
but no interactions. 

On the other hand, lattice methods are not the only available numerical 
approaches for strongly coupled many-body physics. Other potentially exact 
techniques such as Green's function Monte Carlo (GFMC) [2 [HIS], Coupled Cluster 
(CC) dl E] and No-Core SheU Model (NCSM) [S] are in the mainstream, especially in 
modern studies of nuclear structure (see e.g. [a[8l|9l[l0l[ni[l2l[T3l[l4l[l5l[l6lll7l|18l 
[T9l[2^). Still, it should be pointed out that GFMC involves the uncontrolled fixed- 
node approximation, and while this has not prevented practitioners from obtaining 
remarkably accurate results, it is regarded by many as intellectually somewhat 
dissatisfying. Uncontrolled approximations aside, none of these methods scales with 
the size of the problem as favorably as current lattice approaches. In addition, lattice 
methods have robust technical foundations and benefit from a continuous stream of 
attention from the Lattice QCD community. 

Apart from explaining the motivations for using and developing lattice methods, 
and setting the broad context, the purpose of this Introduction is to establish the 
scope for the rest of the article. We have decided to focus our attention on methods 
for non-relativistic systems, more specifically for ground-state calculations. There 
are several excellent books and reviews on relativistic systems, particularly on those 
with gauge degrees of freedom (see e.g. Refs. dU [HI H^), and the topic of finite 
temperature merits a separate article by itself (see e.g. Refs. [231 We shall also 

avoid the all-important issues of time-dependent response functions and their analytic 
continuation to real time, which are topics that deserve their own review articles (see 
e.g. [261 EZ]). 

The motivation for the investigation of lattice methods specifically for non- 
relativistic systems is at least twofold. On one hand, the behavior of the vast majority 
of atomic and condensed-matter systems is governed by electronic Hamiltonians in 
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which relativistic effects enter only as a small correction. The symmetry group is 
thus not given by the Poincare set of Lorentz boosts and spacetime translations, but 
rather by rotations and Galilei transformations, or their discrete counterparts when 
a physical lattice is present. It should be stressed, however, that the richness of low- 
energy condensed-matter physics allows for an immense array of Hamiltonians whose 
characterization defies the above simple-minded binary approach. At this point in 
time perhaps the paradigm of this situation is given by graphene, whose low-energy 
dynamics is controlled by a (relativistic-like) linear dispersion relation plus a (certainly 
non-relativistic) Coulomb interaction. While a few groups are currently using lattice 
field theory methods to elucidate the effects of strong coupling in graphene (see e.g. 
Refs. |28ll2S[30ll3ll[3a|33l|32), we shall not discuss this system beyond this point. 

The other main motivation comes from low-energy nuclear physics, where the 
description is predominantly given by non-relativistic actions, as the mass of the 
nucleons (^ 1 GeV) is much larger than the mass of the low-energy excitations 
(m^ ~ 0.1 GeV), which are pseudo-Goldstone bosons associated with the spontaneous 
breaking of an approximate chiral symmetry (see Ref. |35j). Lattice calculations of 
nuclei are still in their infancy, but are currently being aggressively pursued by the 
Bonn- Raleigh group 1371 [35] , having performed calculations of the lightest nuclei 
up to [jglgn], including a first study of the Hoyle state gTHH]. Ahhough they 
are beyond the scope of this review, the efforts of the NPLQCD collaboration [151 Hi] , 
dedicated to ab initio nuclear physics starting from Lattice QCD, should not be 
overlooked in this context. 

Finally, it should be pointed out that there are dynamic synergistic activities 
among the various areas mentioned above. Indeed, with the advent of ultracold atom 
experiments in the last few years (see Ref. [JSl 1321 HZ] for reviews of the current 
experimental and theoretical situation), the nuclear- and atomic-physics communities 
are working in close collaboration with each other. This applies in particular to the 
study of ultracold Fermi gases close to a Feshbach resonance. The latter are effectively 
spin- 1/2 systems whose interaction can be tuned to vanishing range Tq and large s- 
wave scattering length ag. This seemingly "unnatural" situation is actually very much 
relevant for nuclear structure in general and neutron matter in particular, where the 
ranges are short (rg ~ 2.7 fm) and the scattering lengths are large {ag ~ — 18 fm). 
The stark difference between the nuclear scales (~ fm) and the atomic ones (^A) 
testifies for the property of universality of the underlying dynamics [35] HSl HO] ■ 

Universality is particularly evident in the case of the so-called unitary Fermi gas, 
which is defined as the limit ^ kprQ ^ 1 <C kpa^ — > oo, where kp = (Stt^tt,)^/^ is 
the Fermi wave number and n is the density. By definition, the unitary limit has no 
intrinsic scales other than its density n and temperature T. As a consequence, every 
dimensionful thermodynamic quantity must come as a power of kp, or equivalently 
the Fermi energy Sp = h'^kp/{2M) (where M is the mass of the fermions) times 
a universal function of T/sp (or equivalently /3/Lt, where /3~^ — kgT and /i is the 
chemical potential). In particular, one of the most sought-after universal functions 
corresponds to the energy, which is typically denoted by ^ and often referred to as the 
Bertsch parameter. 

For the same reasons, the unitary Fermi gas displays a number of dynamical 
symmetries and exact relationships. Indeed, as shown in Refs. [5T1|32], the unitary 
limit respects a non-relativistic conformal symmetry (Schodinger algebra [S3J [53]). 
This symmetry is directly related to an exact relation between the scaling dimensions 
of the primary operators of the non-relativistic conformal field theory (which is defined 
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in free space) and the spectrum of the system in a harmonic trap [55l |56l \57\ [58] . 

Aside from those relations, vahd in the hmit of infinite scattering length, there 
are a set of exact universal relations, now widely known as the Tan relations 'SHI HOI 
inn HH ISSl IMl ESI ESI ISZ] (see Ref. [55] for a recent review), that encode the physics at 
short distances and which are due to the short-range nature of the interaction. The 
latter property allows one to write down explicit relations in terms of the probability of 
finding two particles of unequal spin at short distances from each other; this probability 
is commonly referred to as the "contact" . At unitarity, the contact constitutes another 
universal property of quantum many-body physics. The Tan relations are valid at all 
temperatures and couplings, and have also been generalized to spatial dimensions 
other than 3, as well as to bosonic and fermionic systems. The precise calculation of 
the contact in most of these cases is still an open problem as of this writing. 

In this review we shall be only indirectly concerned with the physics of the unitary 
limit, mostly as a test case for the various methods. While this limit is easy to define 
and implement in few- and many-body calculations, it has remained challenging to 
arrive at precise numerical answers for various quantities such as the Bertsch parameter 
and the contact. This is true not only of analytic methods, most of which break down 
in one way or another when approaching unitarity, but also for purely numerical 
methods. For these reasons, this problem is currently one of the most important 
benchmarks for few- and many-body methods in atomic and nuclear physics. 

The rest of the article is organized as follows: in Sec. [2] we discuss the general 
process of formulating the many-body problem on a lattice, and give a first glimpse 
into improvement procedures. In Sec. [3] we present a sample of the current algorithms 
on the market for generating field configurations to be used in the stochastic estimation 
of the path integral. Sec. |4] demonstrates some of the various observables which may 
be calculated on a given set of configurations, focusing on ground-state observables. 
Sec.[5]then discusses difficulties associated with the extraction of the ground-state, and 
provides some techniques for dealing with these problems. In Sec. |6]we consider the 
issue of systematic errors induced by the use of a finite set of lattice points. Finally, 
we present a small selection of results obtained using the methods outlined in this 
article in Sec. [7j followed by our Conclusions (Sec. [s]). An appendix reviewing the 
concept of group theory on the lattice is also provided (App. |9]). 
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2. General considerations 

The main objective of this section is to develop the connection between textbook 
many-body formalism and the expressions commonly used in actual numerical 
calculations. Ultimately, the goal is to develop expressions that allow us to evaluate 
expectation values of quantum operators O in a stochastic fashion, such that, up to 
systematic effects (finite volume, finite lattice spacing, etc.), the exact result is given 
(schematically) by 

{0) = ^J2^nO„ (1) 

n 

where Z = J2n is a positive semi-definite probability measure whose form 

depends on the problem under consideration, though generally not on the observable 
O. Here, n parametrizes the values of a set of auxiliary variables (typically fields that 
live in spacetime), and it is the sum over their (usually unfathomably large) domain 
that is performed using Monte Carlo methods. If one succeeds in doing so, the result 
is then given by 

(^)^^E0„, (2) 

* rs=l 

where Ng is the number of samples of the auxiliary variables, distributed according 
to and the uncertainty in the approximation is of order l/^/Wl if the samples are 
uncorrelated and assuming the central limit theorem holds. The quantity 0„ is the 
expectation value of O in the n-th auxiliary configuration. 

Naturally, limitations arise in the kinds of problems that can be treated in this 
fashion, as it is not always possible (or at least not easy) to define a positive semi- 
definite Vn- This is of course the well-known sign-problem that appears in many 
problems of interest. We shall return to these and related issues below, and restrict 
ourselves to cases where the sign problem is absent. 

Rather than focusing on specific observables, it is often more useful to attempt 
to re-write the partition function Z as a sum of positive terms, as mentioned above. 
As we shall see, in many cases this sets a natural framework for the calculation of any 
observable. 

2.1. Starting point: the partition function & the transfer matrix 

As mentioned in the Introduction, we shall focus on ground-state approaches for non- 
relativistic systems. However, it is convenient to set up the formalism by considering 
the general case of non-zero temperature T. In the grand-canonical ensemble, the 
thermodynamic properties of a quantum many-body system are fully specified by the 
partition function 

Z = e-^^ = Tie-^^^-^'^\ (3) 

where Q. is the grand thermodynamic potential, the trace is over all multi-particle 
states, i.e. over the Fock space, and j3 = {k^T)~^. In second-quantization language 
the Hamiltonian is given, in general, by 

and the particle number is simply 

N = aia^, (5) 
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where sums over repeated indices are iniplied, and aj^ (a^) is the creation (annihilation) 
operator of particles in state a, for a given choice of single-particle states \a). 

In Eq. @ we have written a rather general Hamiltonian including not only the 
conventional one- and two-body terms, but also three- and higher-body interactions. 
We take the point of view that the Hamiltonian in question is effective, in the sense that 
it describes the physics of the problem below a certain scale A, and all the processes 
above that scale are encoded in the so-called low-energy constants (A), and higher- 
body interactions [521 [ZO] • This is expected to be a practicable approach when there 
is a clear separation of scales in the system. While this is the philosophy of modern 
effective field theory (EFT), with few exceptions most Monte Carlo calculations include 
only two-body interactions. Such a simplified theory can be an accurate description 
of dilute atomic systems or neutron matter in the crust of neutron stars (where the 
range of the interaction is small compared to the inter-particle distance), but fails at 
higher densities as well as for nuclear structure calculations. Concerning the latter, 
we expect vigorous research in this direction in the coming years, as the application 
of lattice methods to nuclear physics (with nucleons and pions as the elementary 
degrees of freedom) is currently gathering significant momentum, as mentioned in the 
Introduction. In this review we shall focus on two-body interactions only, but we shall 
comment on the difficulties of implementing higher-body interactions below. 

Following the usual route (see e.g. Refs. [7T1[72]), one may introduce a coherent- 
state representation and go on to the language of fields (with Grassmann fields for 
fermions and complex fields for bosons), such that the partition function of the problem 
at hand becomes 

V^j'fV^ e-^^'f'^ ''''I (6) 

From this point on we shall focus on the fcrmionic case, such that ip denotes in general 
a multicomponent Grassmann field ip^, where a — 1, . . . , iVj, and iV^ is the number 
of fermion flavors present in the problem. The Euclidean action is given by 

S[i;\^]^ J dTd'x{iji{dr^fi)^P^+n[^\^]), (7) 

and the Hamiltonian density reads 

n - vl^o^a + ■ • • , (8) 

where we have implicitly summed over the flavors a, and the interaction terms are in 
the ellipsis. Here, Hq is a differential operator corresponding to the coordinate space 
representation of h^^p in Eq. The full Hamiltonian is of course given by 

[ d^xn[i^'',ij]. (9) 



The fields are defined in all of spacetime, where it should be kept in mind that 
at finite temperature the time direction is compact, extending from r = to r = /3. 
Fermionic (bosonic) fields obey anti-periodic (periodic) boundary conditions at the 
endpoints of the temporal direction. This is the starting point for the formulation of 
the problem in terms of lattice field theory. 

We have thus started from a physical perspective by considering the grand- 
canonical partition function, and moved on to the language of fields, which provides the 
connection with modern methods often seen in Lattice QCD. In spite of the (sometimes 
considerable) technical differences, all of the methods we shall discuss here, and the 
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vast majority of the ground-state methods currently in use, actuaUy share a common 
essential aspect. Indeed, they all rely on the so-called "power" method to extract from 
a given matrix (here the transfer matrix, defined below) the largest eigenvalue and 
eigenvector (here the ground state) (see e.g. Ref. [75]). The idea of starting with the 
finite-temperature formalism and taking the zero-temperature /3 — oo limit is indeed 
nothing other than taking powers of the transfer matrix 

T, = e-'", (10) 

as a filtering process to systematically project out the excited states, assuming that 
one starts with a trial state that has a non-zero (and ideally very large) projection on 
the true ground state. In the language of fields this amounts to the determination of 
the mass gap (or in general, spectral properties) of the fermion propagator by selecting 
an optimal source. 

We shall return to this idea below. In order to appreciate the connection between 
the various methods, which typically come from different areas, it is important to be 
able to jump between the "fields" and "transfer-matrix" languages. We anticipate 
that for many readers this connection will seem trivial. However, it is our experience 
that this is not the norm. 



2.2. Putting the problem on the lattice 

The next step is to discretize spacetime into N!^ x points, and define fields on such 
a lattice. Throughout this review we shall consider a hypercubic lattice with lattice 
spacings bs,br = corresponding to the spatial and temporal directions, respectively, 
as that is by far the most common choice. We denote with tp^ ^ the Grassmann field 
defined on all the space-time points n, r, and tpl^ ^ denotes its hermitian conjugate. 
The lattice Hamiltonian is then 

i7 = ^V4[^^o]n,mV'm + ---, (H) 
n,m 

where n, m are points on the spatial lattice, and [i?o]nm ^ particular choice of a 
discrete representation for the continuum Hamiltonian Hq. In most cases relevant for 
this review, the Hamiltonian will be time-independent. We have again omitted the 
interaction terms, which we turn to below. 

In the absence of interactions, the dispersion relation for a non-relativistic system 
is of course E = /{2M), where M is the fermion mass. The choice of a discretized 
derivative, to represent the momentum operator in coordinate space, determines the 
form of the Hamiltonian of the non-interacting system on the lattice. The various 
choices differ in their approach to the continuum limit, i.e. in their behavior at 
high energies. For the Laplacian operator, the simplest choice is given by the second 
difference formula 

^/J^ E ^[/j+e.+/j-e.-2/j], (12) 
fc=l,2,3 * 

where is the unit vector in the fc-th direction. This yields the correct behavior at 
low energies, as it gives the following kinetic-energy spectrum 

^(p) = E - (13) 

* fc=l,2,3 
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where p ~ hq/b^, and g^. = takes on discrete values as rt^ — —N^/2, . . . ,N^/2. 

Clearly, the above reproduces E{p) ~ ^ for q-^ ^ 1. 

One may reproduce the correct behavior at all energies by using a definition in 
momentum space: 

where 

/n - E ,,3/2 (15) 

is the discrete Fourier transform of /. This yields 

2M 

by definition, at all energies 



E(p) = — , (16) 



The form of Eq. ( 12 ) has the advantage of being as local as possible, in the sense 
that the calculation of the Laplacian at point j requires input only from its nearest- 
neighbors j ± Efc. This means that it can be performed very quickly on a computer. 
It should be kept in mind, however, that the memory layout in current computers is 
linear, such that nearest-neighbors in different directions can be physically very far 
from each other. This fact can affect performance considerably. 



The form of Eq. (14), on the other hand, is as non-local as possible, but has 
the benefit of being closer to the continuum limit than any other implementation. 
References [TH [7S] have analyzed the properties of these choices in the presence of 
interactions using a dynamical mean-field approach (see also Ref. [76 ) 

This is a first glimpse at what is commonly referred to as "improvement" . The 



action resulting from Eq. (14) is an improved version of that of Eq. (12), in the sense 



that it is closer to the continuum limit. In fact, Eq. (14) would lead to what we call 



a "perfect" action for the non-interacting system. We shall return to the issue of 



improved actions in Sec. 2.6 Other forms of improved dispersion relations, which in 



some sense interpolate between cases of Eq. ( 12 ) and Eq. ( 14 1, have been described in 
Ref. [22]. 

The issue of choosing a proper discretization for the temporal derivative dr 
requires some care. Indeed, as noted in Ref. [75], one should be careful when 
transferring fermion operators to the lattice. In particular, the generator of time 
translations, at finite temporal lattice spacing 6^, will not be the discretized version 
of the Hamiltonian shown above; rather, we have 

=-(l-i?,)V,_i, (17) 

where 

= exp(-6,iJ(r)), (18) 

for a general time-dependent Hamiltonian H. Therefore, the generator of translations 
is an effective Hamiltonian given by 

b^H^^ir) = l-e-''^"^-'\ (19) 

Taking this into account, the fcrmionic part of the Euclidean-time lattice action, 
in the absence of interactions, will be 

5 = 5]Vl[ifoU^b (20) 

a,b 
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(we use collective indices a, b to denote points on the spacetime lattice), where 



Ko 



( 1 

-So 




V 





1 








-So 1 

-Sn 



So \ 






(21) 



and where the entries are blocks of size iV^ x A'^, with 1 being the unit matrix 
and = exp{—b^{Ho — fi)). Notice that this discretization amounts to a one-sided 
temporal derivative, which should be contrasted with the relativistic case, where a 
symmetric discretization is used. 

Thus far, aside from the fact that we have only considered fermionic degrees 
of freedom in the problem, the formulation we have described is potentially general 
enough to cover a large number of problems in condensed-matter and nuclear physics. 
Unfortunately, once interactions are taken into account many of those problems are not 
easily tractable with current methods. Special techniques are needed to even begin to 
consider certain cases (e.g. systems with repulsive interactions) due to the appearance 
of the sign or signal-to-noise problem. We shall therefore restrict our attention even 
further whenever necessary, but we will return to signal-to-noise issues below, as it is 
an active topic of research in all areas. 

In order to proceed, we shall focus on two-body interactions, such that in general, 

n,m,j.k 



H = 



V, 



2JnmjkV'j^k 



(22) 



In fact, this problem is already challenging. Indeed, while taking the exponential of a 
one-body operator is relatively easy, it is much harder to do so for a two-body operator. 
Similarly, the only field integral we know how to handle is the Gaussian case. It is for 
these reasons that auxiliary fields are introduced in order to deal with the interaction. 
Unfortunately, even for the simplest interactions there are cases where a sign problem 
arises, such as the case of the point-like repulsive potential. For a purely attractive 
potential, on the other hand, we may "decouple" the two-body interaction by using 
the Hubbard-Stratonovich "trick" , which we now turn to. 



2.3. Hubbard-Stratonovich transformations 

Almost every lattice calculation is based on the fact that two-body interactions (from 
the Hubbard model to QCD) can be represented by mediation through bosonic fields. 
In QED this field is the physical photon field, in QCD it is the gluon field, and 
in theories with contact interactions (Hubbard model, Nambu-Jona Lasinio, Gross- 
Neveu, etc.) one uses a non-propagating "auxiliary field". This auxiliary field is 
introduced using the so-called Hubbard-Stratonovich (HS) transformation (in any 
of its many incarnations) in which the four-fermion operator of the interaction is 
effectively replaced by a path integral over a field (j) that couples to a two-fermion 
operator [751 HOI HI]- Although the HS transformation is an extremely useful step 
toward putting the problem on a computer, it should be pointed out that there is an 
alternative that is particularly powerful for the case of contact interactions, namely 
the worm algorithm (see e.g. Refs. [5^ 15^ 151] ). 
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To fix ideas we sliall consider the case of spin- 1/2 fermions with a contact 
interaction, such that the second term in Eq. ( 22 ) is given by 



(23) 



As mentioned above, we wih focus on time-independent potentials, such that we may 
evaluate the above expression at any time slice r. The HS transformation is based on 
the fact that, at each point a = n, r in spacetime, we may write 

1 



(24) 
(25) 

(26) 
(27) 

where the function represents a contribution whose form depends on the specific 
choice of HS transformation. For these three possibilities (by no means an exhaustive 
list), we have 

Un^'^'^" Continuous, unbounded (Refs. 

Uni^i-'^ + <l>n)^i'^ - K)) Continuous, bounded (Refs. [S3 [77]) 



A discrete version is also available: 

^ 0=±1 

and a continuous but compact form: 

Implementing an HS transformation, the partition function becomes 



Discrete (Ref. [SI]). 



(28) 



The fermionic part of the Euclidean-time lattice action will be 

a,b 

and the dynamics of the fermions is fully specified by the matrix K, which has a 
general structure given by 

/I 

1 

~BM 1 



K\. 



6 





V 



-B 



N -21 





B 



J 



where the entries are blocks of size N!^ x N!^. Notice that the structure of this matrix 
is such that each block column (and row) corresponds to a particular temporal slice. 
It should be clear from the above that the derivative with respect to imaginary time 
that appears in the action is represented here as a one-sided difference. Also clear in 
this equation is the anti-periodic boundary condition along the time direction. In fact. 
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for ground-state studies, the upper-right block Bpf is in principle at our disposal, as 
it merely implements the boundary condition in the time direction, which is infinitely 
long at zero temperature. For instance, Refs. [HSIIHZIIHHI use open boundary conditions, 
i.e. Bj^ = 0; we shall elaborate on this choice in the Algorithms section. 

The form of the i?^ matrices, which in general are functions of 4>, depends on the 
form of the interaction as well as on the choice of HS transformation. We shall see a 
specific example below. 

Since the resulting action is quadratic in the fermion fields, we are now in a 
position to perform the path integration over the fermion fields, which leads to 

J Vc^V[4>], (30) 

where 

V[<f\=p[cf\detK[<f\ (31) 

is the probability of the specific configuration (p. Here and in the following we 
shall neglect overall multiplicative constants in partition functions whenever possible 
(e.g. when integrating out fermions or performing gaussian integrals of any kind). 
(Extensive explanations on gaussian integration over various fields and dimensions 
appear in Ref. [72] ) 

This formulation, which results from the grand-canonical approach to quantum 
statistical mechanics, is the most natural one in connection with finite-temperature 
studies (where the length of the time direction is proportional to the inverse 
temperature) as well as with its relativistic counterpart. 



2.4- The canonical approach 

As particle number is conserved in non-relativistic quantum mechanics, it is often 
natural to start from the canonical ensemble, which leads to a formulation based 
on a Slater determinant trial wavefunction for N particles and leads to a "partition 
function" (technically this is not a partition function as we are not truly summing 
over all possible quantum states, but it is a kind of probability sum) of the form 

Zj, = j V<j>Vj^[<j>], (32) 

where 

Vj^[.p]=p[cf>]detKj^[(f>], (33) 

[K^]^^ = {a\B^^---B^B^B,\l3), (34) 

and ja), a = 1, TV are the single-particle states in the Slater determinant. One may 
define at most linearly independent single-particle states, such that N = is 
the maximum number of particles we may place on the lattice, as befits fermions. In 
general, however, one takes N <^ N!^ to avoid lattice spacing effects (we shall return 
to systematic effects in general in Sec. [6]). 

Ground-state properties are obtained by extrapolation from finite /3 (i.e. finite 
A^^) to the /? — )■ oo limit by fitting to analytic formulas. The performance of 
calculations and subsequent extrapolations in this approach depends strongly on the 
choice of the single-particle basis and trial wavefunctions. For example, Ref. [59] uses 
functions that include some degree of pairing correlations, as does Ref. [UOl, while 
Ref. [5T] uses simple Slater determinants of free-particle states. 
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These two formulations (the canonical and the grand canonical of the previous 
sections) are of course related in more than one way, starting from the fact that the 
matrices that enter the imaginary time propagation have exactly the same form 
given the choice of dispersion relation and HS transformation. In the Algorithms 
section we shall elaborate on the differences and similarities between these two 
approaches. 

As a particular example, if we choose the discrete form of the HS transformation, 

then 

[5fe]„,„' - [^V,]„,„„ (35) 
where the kinetic-energy factor is diagonal in momentum space and given by 

[^]p.p'=e-'^^^^/<'''^'^^^P.p' (36) 
and the potential-energy factor is diagonal in coordinate space and reads 

[Vfc]n,„' = (l + ^'/'(fc,n)),5„^„,. (37) 

where C = e^^^ — 1. (Notice that if g is negative, which corresponds to a repulsive 
interaction, then C is negative as well, and its square root is complex. This brings 
about the sign problem, as mentioned before.) 

As the path integral is to be performed over (f>{k,ri), the actual form of the 
operator V^, will obviously vary wildly throughout a calculation. The operator T, on 
the other hand, is constant. Therefore, in practice, T is computed only once at the 
beginning of a calculation and stored in memory. Furthermore, it is best to compute 
it in momentum space where it is diagonal (it is a fully dense matrix in coordinate 
space) and apply it as needed by using fast Fourier transforms (FFT). 

The strategy of using FFT to apply operators that are non-local in coordinate 
space but diagonal in momentum space, is commonly referred to as "Fourier 
acceleration" (see e.g. [92l [93l |94l [95]). To the uninitiated this may appear to be 
a ridiculously convoluted way to proceed. Why not do everything in coordinate space 
instead of switching back and forth? The reason is that applying a dense matrix of 
dimension M x M to a vector is an operation whose computational cost scales as 
M^, whereas using FFT for this purpose changes the scaling to MlogM. Moreover, 
modern FFT libraries, such as FFTW |96j are so advanced that their performance 
is very difficult to improve, even in cases where the matrix in question is somewhat 
sparse in its coordinate-space representation (e.g. when using a nearest-neighbor finite- 
difference formula for the momentum operator). 



2.5. A comment on many-body forces 



As mentioned above, in this article we restrict ourselves to two-body forces. Higher- 
body forces, however, play a crucial role in nuclear physics and should therefore be 
somehow accounted for. In Ref. [33] a strategy is used to include interactions beyond 
two-body in a perturbative fashion. 

If the perturbative approach is justified, i.e. if the couplings (assume n > 2) 
to the interaction operators JT„ are indeed much smaller for the higher-body forces 
in question than for the two-body sector, then one may expand the action around 
vanishing g^, such that the partition function becomes 



lg„=0,ri>2 " 



n>2 



dZ 
dgn 



n,m>2 



dgndg„ 



.(38) 
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In practice this is performed introducing sources coupled to one-fermion operators. 
The evaluation of observables will therefore include propagators reflecting the insertion 
of the interaction. A different way to proceed is to attempt to perform an 
HS transformation on the many-body interaction in question. Unfortunately this 
introduces a sign problem in almost any case of interest. 



2.6. Improved Actions & Transfer Matrices 

At this point we have presented the various steps necessary to connect the formalism 
of finite-temperature many-body quantum mechanics with the expressions amenable 
to numerical calculation. Before moving on to specific algorithms to sample the 
probability measure, and then to the problem of computing observables, it is 
worth pointing out a recent development which modifies the zero-range Hamiltonian 
presented above in order to allow faster approach the continuum limit. This represents 
a kind of improved action. Improved observables can also be defined using the same 
methods, as we shall see in a later section. 

The starting point is the transfer matrix T. By definition, as we saw above, 

r = exp{-b,H) ~ / d(t,{k)lCVM + 0{bl), (39) 



where /C and Vj, are as in Eqs. (361 and (37), respectively, and we have set /x = 



In order to reduce lattice-spacing effects in non-relativistic systems Refs. [^3 US] 



employed the following strategy. One may promote the constant C in Eq. (37) to 
an operator C{p) that is diagonal in momentum space, such that a more general 
HS transformation is defined. In order to fix the form of the function C{p), one 
parametrizes its form using a convenient set of basis functions 02n{p), 

. No-l 

^(P) = ]i^ E C'2„02„(P), (40) 

n=0 

and then diagonalizes the transfer matrix in the two-particle subspace of the Fock 
space, which is given in the momentum representation by 



where ^ = A'^^ is the lattice volume and T{k) = (fc^ + fc|)/(2M). 

The exact eigenvalues of the transfer matrix corresponding to a given choice of 
scattering parameters is given by Liischer's formula [Ml f99j : 

pcotSM = + lr,p' + 0{p^) = \s{r^) (43) 
where a is the scattering length, r is the effective range, 77 = and 

5 (77) = lim y ~ f ^ - 47rA, (44) 

n ' 

where the sum is over all 3D integer vectors. 

In the case of the unitary limit, where by definition 

pcot<5o(p) = 0, (45) 
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Figure 1. (Color online) Plot of pcotS(p) as a function of r}'^ = iJj/Po (where 
Pq = 2it/L), for Nx = 20, = 0.05, and levels of improvement A''^ = 1 — 5, for 
a Galilean invariant definition of the transfer matrix. 



the above procedure yields considerable improvement in the approach to the 
continuum. Indeed, with a single coefficient Cq one may tune the scattering length 
a, but the effective range r remains finite due to lattice-spacing artifacts (systematic 
effects will be discussed in more detail in Sec. |6]). Naturally, by introducing more 
parameters one may tune the effective-range expansion in Eq. ( 43 ) to higher accuracy 
(see Fig. [T]) , thus systematically eliminating the need for extrapolations to the dilute 
limit and leaving only finite- volume effects unaccounted for. 

Thus, tuning to the unitary point requires knowledge of the roots of 5'(ry). The 
first few roots are quoted in the Appendix. The result of the fit is shown in Table [TJ 
for various improvement levels. 
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Table 1. Results of fitting the coefficients C„ to the low-energy spectrum of the 
two-body problem at resonance, in a box of side = 16, for an imaginary time 
step = 0.05, in lattice units. 

Cq Ci C2 

~i 0.68419 = = = = 

2 0.53153 0.07896 _ _ _ 

3 0.49278 0.04366 0.01807 

4 0.47217 0.03711 0.00784 0.00467 

5 0.45853 0.03331 0.00718 0.00132 0.00129 
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3. Algorithms 

In the previous section we outlined the general strategy that is typically followed to 
write down a two-body Hamiltonian in a form that can be treated with Monte Carlo 
methods. Our main development was expressing the partition function as a sum over 
auxiliary-field configurations: 



where 'P[(/>] generally involves the determinant of a large matrix whenever fermions 
are involved. 

Once a suitable probability measure Vlcj)] has been chosen, one faces the task 
of constructing an efficient method to sample uncorrelated field configurations that 
obey such a probability. Here, the importance of the word "efficient" should not be 
underestimated. Indeed, the possibility of evaluating a path integral stochastically is 
a remarkable step forward given the size of the domain of integration, but stochastic 
integration of quantum systems would not be viable without methods for fast sampling 
of V[(j)]. This consideration is quite general, applying to relativistic as well as non- 
relativistic problems. Indeed, none of the endeavors of the Lattice QCD community 
would be possible without fast sampling methods. In this section we describe some 
of these methods, specifically for the case of non-relativistic systems. The problem 
at hand is equivalent to that of constructing an efficient random number generator, 
where the output is not just one number but a field defined in space-time. 

Among the methods available to tackle this problem, one may identify two 
large classes: the so-called heat-bath methods, which allow for the direct generation 
of uncorrelated samples; and Markov-chain methods, based on the Metropolis 
algorithm |100j . While the former are extremely useful when studying for instance 
pure gauge theories [1011 1102L 11031 1104) , their use is much more limited in the case 
of theories that involve dynamical fermions due to the complexity (non- linearity and 
non-locality) of the fermion determinant present in V[(j)]. As we shall see, however, 
heat-bath approaches are still useful for non-relativistic fermion theories in certain 
formulations. Markov-chain approaches, on the other hand, are the most popular and 
generally viable when dealing with fermions (as long as there is no sign problem), and 
we shall therefore discuss them first. 

3.1. Determinantal Monte Carlo 

One the simplest methods to sample a non-trivial fermionic probability is known by 
the name of Determinantal Monte Carlo (DMC). In this algorithm, one generates a 
sequence of samples of the auxiliary field (j) (the Markov chain) by performing localized 
changes, typically of random magnitude affecting a random set of locations. The 
new configuration, (/>new' is then accepted or rejected according to the Metropolis 
algorithm, i.e. with a probability of acceptance given by 



where V[(l)\ = detK[(j)\. 

As successive configurations are strongly correlated with each other, snapshots of 
the field should be taken at sufficiently long intervals along the chain. Such a process 
of decorrelation can be time consuming, depending on the proximity to critical points 




(46) 




(47) 
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[]. The minimal requirements for this algorithm to function and yield the correct 
probability measure are ergodicity and detailed balance (though technically a weaker 
condition is enough, see e.g. Ref. |105| for a nice overview). Ergodicity implies 
that long enough chains will eventually touch every point in the space of possible 
configurations. Detailed balance, i.e. ensuring that a given update has the same 
probability as the reverse process, implies that the Markov chain has as a fixed 
point. 

When moving from one configuration to the next, it is possible in some cases 
to update clusters of points at a time, as explored many years ago in the context of 
the Ising model |106l 11071 1108j . However, such attempts to take long random steps 
in configuration space in our case will in general, though not always, result in a low 
probability of acceptance. One is therefore forced to take relatively small steps. On 
the other hand, taking very small steps should be avoided as well, as they yield high 
acceptance but slow decorrelation and are therefore very ineflticient. The rule of thumb 
is then to compromise by choosing updates such that the acceptance rate is in the 
range 40-60%. 

The DMC method is not only among the simplest, but is also one of the oldest and 
remains in wide use (see e.g. |109j for a recent use of DMC); it is also very effective, but 
unfortunately not particularly efficient. The main drawback of the "plain" version of 
DMC presented here is that localized updates of the field represent a huge cost from 
the computational point of view. Indeed, strongly correlated systems are generally 
not amenable to global modifications of the field, and the (integrated) size of the 
randomly updated regions generally does not scale up with the volume of the system. 
As a consequence, the cost of a full update (a sweep), which is a measure of the 
cost of obtaining a useful configuration, scales as NN^V^, where N is the number 
of particles, is the number of lattice points in the imaginary time direction, and 
V = is the spatial volume. This estimate results from the evolution of N particles 
in imaginary time over A^^ steps, with a cost per step equal to that of dense matrix- 
vector multiplication, which on a basis of size V scales as V^. The cost of local 
updates contributes an extra factor of N^V, because the auxiliary field lives in all 
of spacetime, which yields the final result. This estimate does not account for the 
possibility of Fourier acceleration, in which case the cost of a time step is proportional 
to V log V, resulting in a cost per sweep of NN^V^ log V. 

3.2. Hybrid Monte Carlo 

As of this writing, the use of Hybrid Monte Carlo (HMC) in non-relativistic physics 
is surprisingly scarce, both at zero as well as finite temperature. While this powerful 
algorithm was developed many years ago in the field of Lattice QCD |110l llllj , its 
application to strongly correlated systems in other fields, in particular condensed 
matter and low-energy nuclear physics is still in its infancy (with some exceptions, see 
Ref. [771 HH])- 

The HMC method, as DMC, is based on the generation of a Markov sequence of 
configurations that are accepted or rejected as dictated by the Metropolis algorithm. 
The fundamental difference between DMC and HMC (in its various forms) is that the 
latter enables global updates of the field, which reduces the cost per sweep to NN^V'^, 
or NN ^V\ogV when Fourier acceleration is possible. This dramatic improvement is 
accomplished by introducing molecular dynamics (MD) into the updating procedure 
in the following way. 
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Given the probability measure V[<j)] and corresponding "effective" action 

S,s[^]^-logV[cl)], (48) 

one introduces a fictitious gaussian-distributed momentum field tt, which modifies the 
partition function only by a multiplicative constant, i.e. in a way that is immaterial 
to the dynamics of the system: 



where 



V^V[<j)]^ J V(t)VTTr[(l3,T:l (49) 
V[<P,n] = exp ( - ^ % = exp{~nMJ^), (50) 



2 



and 



^MD^E^ + ^cffl-^]- (51) 



n.r 



Since we have not modified the true dynamics of the problem, this new probability 
P[(j),Tr] is physically equivalent to the original one. The reason it is advantageous to 
work with this enlarged problem is that one may use MD as the updating strategy to 
explore configuration space. Taking a gaussian distributed tt and any choice of (f> as 
the initial conditions, one solves the classical evolution equations dictated by ^j^p- 



<, = f^„,,[0]^-^|«M, (52) 

where the dots indicate a derivative with respect to the fictitious MD time ijyjQ • Notice 
that these equations are global, i.e. they affect each and every point in space-time, 
such that upon integration along a trajectory of length Ty^j-^ 1 one obtains a fully 
updated field configuration, i.e. one has performed a sweep. If these equations were 
integrated exactly, the fictitious energy given by 'Hy^Y) would be conserved, and the 
probability of acceptance of any trajectory would be exactly 1. It is this combination 
of MD with the Metropolis algorithm that gives the "Hybrid" in HMC. In practice, 
integrators are precise but not exact, such that a Metropolis step, comparing the 
probabilities at the beginning and at the end of the trajectory, is needed to guarantee 
that the correct probability distribution is being sampled. Resampling of the fictitious 
momentum tt at the beginning of each trajectory is needed to actually perform the 



integral over tt that appears in Eq. ( 49 ) . This also takes the random walk in a new 
direction and into a different MD classical orbit, which favors faster decorrelation. 
In order to maintain detailed balance, as required by the Metropolis algorithm. 



the MD equations (52) must be integrated in a reversible and area-preserving fashion 
(see e.g. Refs. [1101 lllll 1105] ). This is accomplished by implementing a symplectic 
and symmetric integrator, such as the leap-frog algorithm, or more sophisticated 
approaches like the Omelyan integrators |112j . 

It should be stressed that, in contrast to the DMC algorithm, the current form of 
HMC applies only to the case in which the field 4> is continuous. This is a prerequisite if 
classical differential equations are to be used for the MD evolution, and it is somewhat 
unfortunate, as there are a variety of formulations that rely on discrete fields, which 
have the obvious advantage of requiring significantly less storage space than their 
continuous counterparts. 
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3.3. Pseudofermions and connection to the finite temperature case 

One of the most common versions of the HMC algorithm, used routinely in Lattice 
QCD, utilizes the so-called pseudofermion fields x to represent the determinant via 
the following relation: 

det Af [0] cx J Vx^Vx exp {-x^ M-\4>]x) , (53) 

where we assume that all the indices are contracted in the exponent on the right- 
hand side, and it is assumed that the dynamics encoded in M is such that the 
gaussian integral is well defined. In practice, this integral is evaluated stochastically, 
exploiting its gaussian form. The name "pseudofermion" arises from the fact that 
these are bosonic complex- valued fields whose dynamics is given by M~^, which 
satisfies antiperiodic boundary conditions in the time direction. To our knowledge, 
pseudofermions have never been used in non-relativistic ground-state calculations, 
possibly because in that case the determinants are substantially smaller than in the 
finite-temperature and/or relativistic cases, such that an exact calculation of det M is 
feasible and even convenient in comparison. 
In this case the MD force takes the form 

F„,,[0] = ™ [M-'m , (54) 

where again all the indices are implicitly contracted. It is easy to see that a fermion sign 
problem will ensue unless M is positive semidefinite. The latter happens, in particular, 
in cases with an even number of flavors, where (under certain circumstances, e.g. spin 
symmetry) we may define a matrix K such that det M = det K . 

A less common version of the HMC algorithm combines the advantages of both 
the DMC and HMC algorithms, without using pseudofermions, resulting in the so- 
called Determinantal Hybrid Monte Carlo algorithm. The idea behind this approach 
is to make use of the fact that 

detAf = det (1 -I- [(?;.] )^/, (55) 

where 

U[(j3\ =Bj^^--- B^B^B^ (56) 

(see Sec. [2]). 



Notice that the determinant on the left-hand side of Eq. ( 55 ) is over a spacetime 
matrix, whereas the one on the right-hand side is purely spatial. The fact that 
this simplification is possible is directly related to the absence of backward temporal 
propagation in non-relativistic theories. 

The zero-temperature version of this discussion is directly connected to the use 
of a Slater determinant jao) as the starting point or "guess" for the ground-state 
wavefunction, which we mentioned in Section [2] In that case, 

V(f>det{U[(j>])^f, (57) 

where the determinant is to be taken over the space generated by the single-particle 
(s.p.) orbitals (pj that make up |ao), such that 

-1 m 



F^M=Nf tr 



(58) 



where, as with the determinant, both the trace and the inverse are to be interpreted 
as restricted to the subspace generated by the s.p. orbitals (pj. 
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3.4- Heat-bath approach 



An interesting recent development in the field of lattice fermion calculations concerns 
exploiting the fact that the non-relativistic ground-state problem can be formulated 
in a way that the fermion determinant is identically equal to 1, for all the field 
configurations [86l Ull UHl HH Ull I113j . Indeed, as mentioned in the Introduction, 
ground-state calculations give us the freedom to choose the form of the boundary 
condition in the time direction, and in particular we may choose an open boundary 
condition, such that the dynamics is governed by the triangular matrix 

f 1 ... \ 

-B^ 1 ... 

1 ... 



Kcb] = 



\ 



-B 



1 



-B 



N-l 





1 / 



(59) 



Not only is det if = 1, but the triangular form also gives us direct access to the inverse: 

/ 1 ... \ 



B2B1 



1 



B, 



(60) 



This constitutes a case in which, trivially, one can sample uncorrelated field 
configurations directly according to the relevant probability measure, and it thus falls 
into the category of heat-bath approaches. This is a dramatic simplification of the 
original problem that comes with advantages and disadvantages. 

On the positive side, a constant probability measure means that the cost of 
generating configurations is essentially insensitive to the number of particles and 
the size of the lattice, depending only on the speed of the uniform random number 
generator of choice. Field generation in this approach is therefore instantaneous 
compared to the other methods. Furthermore, a single set of configurations and 
propagators may be generated and then used for any number of observables desired. In 
addition, a constant probability implies no sign problem, regardless of spin asymmetry, 
nature of the interaction (attractive/repulsive, two-body, three-body, etc.), etc. This 
approach therefore constitutes an intriguing possibility for ground-state studies of a 
variety of systems in condensed-matter and nuclear physics, where the sign problem 
has been a major roadblock to progress. 

On the negative side, this formulation is restricted to zero temperature, where one 
may legitimately ignore the (anti-periodic) boundary condition in the time direction. 
Furthermore, a constant probability implies that importance sampling is not possible: 
all the configurations are in principle equally important. In practice, the sign 
problem reappears here in the form of a statistical "overlap problem" : the random 
configurations may not have much weight for the observable of interest. Using a very 
large number of configurations (^ 10*) and a sophisticated statistical analysis (see 
Sec.[5| it is still possible to perform calculations and obtain results with unprecedented 
precision. 

Perhaps one of the main advantages of this method is not that it "solves" the 
sign problem, but rather that it replaces it with a new problem that we have a better 
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chance to tackle, in particular if aided by physical insight into the system at hand. 

3.5. Open-ended imaginary-time evolution, re-weighting and branching random walks 

Extracting ground-state properties entails performing calculations for j3 as large as 
possible. We will return to the issue of extrapolating to the /3 — >■ oo limit from finite 
(3 data below. It is important to underline, however, that serious technical difficulties 
arise at large f3. As explained most recently in Ref. [97l I113j . and as we review 
below, an interval of exponential approach to the asymptotic ground-state properties 
(even for simple observables such as the energy) is typically followed by a dramatic 
deterioration in the signal-to-noise ratio, which may render the whole extrapolation 
effort dubious, or at least quite complicated. 

The algorithms explained in the previous subsections generally involve fixing the 
spatial and temporal extents of the lattice, followed by sampling configurations of the 
auxiliary field defined on that lattice. Notice that the heat-bath algorithm, however, 
does not require a temporal direction of fixed extent, as the probability measure does 
not depend on the size of the lattice in any way. 

There is a different set of algorithms that have been in use since at least the mid- 
1990's, which allow for an open-ended temporal direction, such that in practice it has 
been possible to reach temporal extents a factor of at least 2 or 3 larger than with the 
methods outlined above. Clearly, following this strategy is at odds with the common 
practice of saving auxiliary field configurations of fixed size. These are, however, 
powerful algorithms that have largely remained unknown in various communities, and 
so we wish to outline them here. 

The starting point is, as above, the factorization and Hubbard-Stratonovich 
transformation of the transfer matrix: 



This is valid for each temporal slice, and we may therefore take a "guess" 
wavefunction |\1/q) and apply a sequence of T's for a given spacetime field configuration 
(/>. The problem remains, of course, of summing over all possible </> fields in order to 
recover the physics we are interested in. Should we choose to pick configurations 
purely randomly, we would obtain the heat-bath algorithm of the previous section. As 
mentioned before, this may lead to substantial statistical noise. 

The idea we wish to outline here, which first appeared in Refs. |1141 I115j . is 
essentially what in certain areas is called "re- weighting" , but with a crucial twist. 
Instead of using a flat probability measure, one may assign a weight to the auxiliary 
field configurations; the crucial point is in introducing this weight by defining a 
probability in which the temporal dependence factorizes completely. This property 
of factorization is trivially true for the heat-bath method, but is certainly not satisfied 
for the DMC, DHMC or HMC approaches, where all of spacetime is inextricably linked 
into the fermion determinant. In this fashion, one may perform an open-ended random 
walk in imaginary time. 

The first step is to consider not a single Slater-determinant wavefunction but a 
sum of them, all with equal probability, making up the initial guess for the ground 




(61) 




(62) 
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state l^'o): 



N. 



^o)=y\vl). (63) 



k=l 



The imaginary-time evolution from step i to i + 1 is then carried out on each 
Slater determinant by applying the HS-transformed transfer matrix: 

= J i?0,+ir[0,+i]|^l). (64) 

In principle, in the above path integral, all the field configurations at that time slice, 
4>t+ij E^re equally important. To improve the efficiency of the temporal evolution we 
introduce the re-weighting strategy mentioned above. 

For each evolved Slater determinant state (what is often called a "walker") \(pl,) 
we define an importance function 

OA^i) = {^tH)- (65) 
With this weight we can now write the imaginary-time evolution as 

= J (66) 

where we have defined a modified walker 

= OA^iM) (67) 

and a probability 

which depends on the initial and final states Iv?*-*-*) and \(f^^~^^'>), respectively. 

At this stage, the algorithm appears to be reduced to sampling the auxiliary field 
4>t+i with the probability 7^/j['/'t+i], and thus construct the evolution operator for the 
modified walkers. It can be shown, as in the cases we discussed before, that in the 
presence of an equal number of fermions for each flavor (and for an even number 
of flavors) this probability is positive semideflnite, and this holds in particular for a 
particle-projected BCS state ^QCT. At the end of the problem the unmodifled walkers 
are recovered simply by dividing by the relevant Oj,. 

There is however a problem in that 7^^. is not normalized, and typically importance 
sampling will introduce a normalization automatically (which is another way of saying 
that a normalized probability is required). We are therefore forced to consider the 
expression 



where 



(70) 



and 



Ni = j (71) 

is of course the normalization of Vi. 



kl 
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The evolution in imaginary time then proceeds by obtaining samples of (pt^i 
that obey the normalized probability measure Vf., and using the stochastic estimate 
of the path integral as the evolution operator: 

P^t+in['^t+i]r[<^t+i] ^ 1^ ET[0?+i]. (72) 

'^'^ 9=1 

The normalization constants corresponding to each walker are separately 
accumulated in the so-called "weight" variables according to 

(73) 

The calculation of the normalization is in general not trivial. However, as 
explained in Refs. [1141 IllSj . this can be overcome by updating one lattice site at 
a time. 

After a period of equilibration, i.e. for large enough t, one can expect to obtain 
an estimate of the ground state, which is given by 

|vl/„)^^^*|^*), (74) 

k 

up to an overall normalization factor. 

In practice, the weights of certain walkers will quickly grow after a few time steps, 
while others will decrease to zero. It is for this reason that the concept of "branching" 
is introduced. The idea is simple: every few steps in imaginary time, consider the 
possibility of randomly keeping or discarding some of the walkers. The number of 
steps (or equivalently the size of the time step) should be adjusted in such a way 
that roughly 40 — 60% of the walkers are discarded at each instance. Of the walkers 
that remain, some of them are used to generate two or more new walkers, adjusting 
the number with a given algorithm, in such a way that the total number of walkers 
remains approximately constant throughout the calculation. At that point it is also 
important to re-orthonormalize the whole set of walkers to avoid stability issues. By 
allowing for the possibility of branching, one is essentially re-starting the calculation 
with a many-body wavefunction that is closer to the ground state than the original 
guess. In this way it is possible to evolve in imaginary time for much longer than with 
other algorithms. 
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4. Observables 

4.I. Calculating the energy 

The simplest observable to compute using lattice methods is the ground-state energy 
of a system. One method for doing so is to construct the A^-body correlation function, 

CnW) = ^J 2?Wtp^e-^['^''^-^lvi/(^)(/3)vI/]^^)(0) , (75) 

where Z is the partition function, and the source (sink) is given by 

f dxi---dXnA^^\xi---Xn)Mxi)---MXn), (76) 



where A^'^^ is an A^-body wavefunction and the fields -04 are chosen with the 
appropriate quantum numbers such that ^'^"^ has nonzero overlap with the state 
of interest. The simplest form for is a product of non- interacting single-particle 
states, however, more sophisticated choices, for example including pairing correlations, 
may be incorporated into the wavefunction to have superior overlap with the 
ground state. This will be discussed in more detail in Sec. [5] 
Upon integrating out the tp fields we have 

Cn{P) = ^ ! V(t>Vm{<l),P), (77) 



where 

C(0,/3) = /(^i(0,/3),52(0,/3),---) (78) 

where / is a function of single- or multi-particle propagators, depending on the form 
of VP. This function should be properly (anti)symmetrized for (fermionic) bosonic 
fields. Using the simplest example of sources and sinks composed of single-particle 
states |q) , |/3), / becomes a Slater determinant of single-particle propagators, so that 
C = AeiSa,i3, where 

= (alBAT^-.-SgSaSil/?), (79) 

and w e rep roduce the canonical partition function, Eq. |32[ The form of the propagators 
in Eq. 79 may be recognized as a product of transfer matrices, = e"*"^^, projected 
onto the chosen initial and final states. 

One may insert a complete set of energy eigenstates into the correlation function. 



Eq. 75 to show that for a given Euclidean time /3, 

m,n r — 1 

= Zoe-^^o -f Zie-'^^i + • • • , (80) 

where Za = {"^ src\o) snk) is the overlap of the a-th energy eigenstate of the system 
with the chosen sources/sinks and Ea is the energy associated with that state. To 
extract the ground state, it is often useful to construct what is commonly called the 
effective mass, 

For sufficiently large /3 ^ E'l — i?o, the ground state energy will dominate as excited 
state contributions are exponentially suppressed. By plotting this quantity as a 
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function of /3, one looks for a plateau to indicate the elimination of excited states. 
The plateau region may then be fit to extract the ground state energy as well as a 
statistical uncertainty using standard methods. 

One may also choose to use the thermodynamic definition of the energy, E(/3) = 
— . For the large-/? limit one obtains an expression similar to that derived for 

the correlator, 



where 



EW) 



bE = 



Taking for example 



0- 



Zo 
the 



En 



5b e 



{El — Eo). 



(82) 



(83) 



canonical partition function obtained using a Slater 



determinant state Eq. (32), the observable calculated on each field configuration will 
be 



-Tr 



(9/3 



(84) 



In Sec. |4.3| techniques for improving observables will be discussed, using this form of 
the energy as a first example. 



As seen in Eqs. (81 ),(82 1, information about excited state energies is also available 
from these methods. Fortunately, there exists a vast literature on sophisticated 
techniques developed by the lattice QCD community which one may borrow to tackle 
the problem of extracting excited state energies [1111 [TTZl [111111111123 [HI]- Recently, 
the low-lying excited states of ^^C have been calculated using non-relativistic lattice 
methods, including a signature for the Hoyle state [ITj . 



4-2. Scattering parameters 

Once the energies of a system have been calculated for a given set of parameters there 
are many quantities of interest which may then be derived. Scattering parameters 
form one class of quantities which can be calculated directly from the energy via 
the Liischer formula, Eq. (43). According to the Maiani- Testa theorem 



Euclidean Green's functions cannot give information about infinite volume Minkowski 
scattering matrix elements except at kinematic thresholds. The Liischer formula gives 
an indirect method for calculating these parameters by relating the volume dependence 
of multi-particle energies to the corresponding infinite volume scattering phase shift. 
As shown in Fig.[2j the energy eigenvalues for two particles in a periodic box are given 
by the intercepts of the corresponding irLp cot S with the function S{ri), defined in 
Eq. (|44]). 

In [124^ it was emphasized that the Liischer formula is valid for any value of the 
scattering length, regardless of the size of the box, so that moderately-sized lattices 
may be used even for scattering processes with unnaturally large scattering lengths. 
The authors of this work derived approximate expressions giving the explicit volume 
dependence of the energies of two-body, s-wave scattering states in the limits of small 
and large scattering length, a/L. In the former limit, the lowest-lying energy is given 

by 



En 



Ana 



1 



(85) 
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Figure 2. •S'(r)) (solid red) and -ivLp cot 5 (dashed) vs. rj'^ = (^27rj ' "^^^ 
TzLp cot S correspond to the scattering parameters ro/a = — 0. 1 , with the following 
volumes: L/\a\ = 2 (blue), L/\a\ = 4 (magenta), L/\a\ = 8 (gold), L/\a\ = 10 
(green). The energy eigenstates for the corresponding volumes are given by the 
intercepts of the curves. 



while for the latter one may use 
47r2 



En = 



di + d2Lpcot So 



(86) 



where ci , C2 , di , ^2 are numerical constants determined in |124j , along with analogous 
expressions for bound states and first excited states. However, in addition to the 
universal power law dependence of these energies on the size of the box, there exist 
non- universal exponential corrections due to finite range effects jl251 11261 , setting a 
limit on how small the lattice may be for a given scattering process at the desired 
accuracy. 



Generalizations of Eq. ( 85 ) to many-body systems have been accomplished using 
perturbation theory for systems of bosonic particles |1271 11281 1129) , and successfully 
applied in lattice QCD calculations of multi-meson states |130U131lll32l[TB51ll34j . The 
Liischer formula may also be generalized in many other contexts, for example, relations 
have been derived for asymmetric boxes |135[ll36j . higher partial waves |137|I138|I139) . 
moving frames |140|, I141j , moving bound states [1421 1143j , and multi-channel processes 

[mills]. 

All of these versions of Liischer's formula begin with the assumption that the 
particles exist within a periodic box, a common situation for lattice calculations. In 
some cases it may be advantageous to consider a different form for the IR cutoff. 
For example, an adaptation known as the spherical wall method can be useful for 
determining phase shifts at higher orbital angular momenta, as well as spin-orbit 
coupling and partial- wave mixing [146] . Here, a hard spherical boundary is imposed in 
position space, and the eigenvalue problem then becomes one which is well-known from 
standard nuclear theory texts: the Schrodinger equation is solved given the condition 
that the wavefunctions must vanish at this spherical boundary. The hard cutoff 
effectively removes the complications induced by multiple-scattering effects arising 
from the periodicity of the lattice, however, additional complications arise due to the 
lack of spherical symmetry on the lattice. This method has proven useful for setting 
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unknown operator coefficients in lattice chiral EFT methods [37l EHl [39l HQ] . 

Another possibihty is to use the eigenvalue solution for two particles in a harmonic 
potential [1471 11481 I149j . Due to progress in the development of nuclear EFT 
approaches using an oscillator basis [1501 11511 1152) and for systems confined in a 
harmonic potential |153U1541ll55j . it has been suggested f |156lll57] ) that the spectrum 
of complex nuclear systems in a harmonic potential may be useful in extracting 
scattering information using the relation 



where e = E/uj. This technique has been demonstrated for two fermions at unitarity 
|158j . and the two nucleon system |157j in the continuum. 

It was suggested in [88] that the use of a harmonic potential in lattice calculations 
may offer some advantages, such as a reduction in noise and fast convergence to the 
ground state due to the confinement provided by the potential. On the other hand, the 
potential introduces new physical scales which may induce systematic errors that must 
be accounted for. For example, because the harmonic potential is non-zero everywhere 
(except for the origin) there will be modifications to the short range behavior. These 
finite range effects will not be exponentially suppressed as they are using the standard 
Liischer relation, so that an extrapolation to free space (w — > 0) may be required. 
Systematic errors associated with harmonic potentials will be discussed in Sec.[6| 

4-3. Improved Ohservahles 

4-3.1. Energy. The work of Refs. [571 ES] tackles the problem of reducing UV 
lattice effects by defining an improved the transfer matrix, using a generalized HS 
transformation and a low-momentum expansion. As we saw at the end of Sec. (2] one 
may then tune the coefficients C„ such that the two-particle energy spectrum matches 
the one required by Liischer's formula for the lowest Nq eigenvalues, given the desired 
values of the scattering parameters. 

The above procedure represents a significant step forward in mitigating lattice- 
spacing effects in MC calculations, especially considering that it requires only a small 
coding investment for its implementation in extant MC codes, and it results in minimal 
computational overhead. 

One may use the same technique to eliminate lattice-spacing effects from 
observables. This is particularly useful in connection with improving finite 
temperature lattice calculations, such as those of Refs. |109| [25l Il59| I160j . Indeed, in 
those calculations, as well as in similar ground-state approaches, the transfer matrix 
is not the only object carrying lattice-spacing effects: the operators used to compute 
expectation values also suffer from the same problems. 

The main idea behind the method remains the same as in Ref. [97) . Liischer's 
formula provides us with the exact two-particle spectrum, such that the eigenvalues 
of the exact two-particle transfer matrix and its derivative(s) are known: 



where E2 are the exact two-particle energies in a continuous box [J] 

X We use 7^ here for didactical purposes, but one may use the symmetric expression 7^^ 7^ instead, 
with straightforward modifications, as long as this same expression is used in the actual Monte Carlo 
calculations for a double step in the imaginary-time direction. 




(87) 



E2exp{-b.^E2) 



(88) 
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Table 2. Results of fitting the coefiicients D„ to the low-energy spectrum of 
the two-body problem at resonance, in a box of side A^^ = 16, for an imaginary 
time step = 0.05, in lattice units. 



No 


Do 


Di 


D2 


^3 




1 


-14.76869 










2 


-11.54894 


-1.74519 








3 


-10.74506 


-0.96946 


-0.40164 






4 


-10.31974 


-0.82605 


-0.17494 


-0.10404 




5 


-10.03874 


-0.74266 


-0.16064 


0.02948 


-0.02878 



In order to match this spectrum, we take the derivative of Eq. (41 1 with respect 
to b^, evaluated between eigenstates of the proposed transfer matrix 7^. We thus 
obtain, using the Feynman-Hellmann theorem, 

d{E\%\E) 



where 



and 



Pr 



= {E\e 



2M 2M 



1 dA[p,) 
2V db^ 



[K, + U,] e- 



MPr) 



E) 



2V 



1 

2V 



E 

n=0 



(89) 



(90) 



(91) 



The rest of the recipe consists in taking the right-hand side of Eq. (89) and 
fitting the coefficients D„ such that the first eigenvalues of the exact expression 
Eq. (88) (which correspond to the lowest eigenvalues prescribed by Liischer's formula) 
are reproduced. We may assume at this point that the coefhcients C„ are known, such 
that the eigenvectors \E) are fixed when we set out to find the In that case, the 



-D„ are determined by a linear system of equations of order Nq 



E 

n=0 



where 



and 



(92) 



(93) 



Ye = Eexp{~b^E) - {E\e'^ K^e^^lE) (94) 

Once the coefficients D„ have been determined, one can use this in a lattice 
calculation simply by replacing A ^(p), and taking dC^/db^ — _D„. The results of 
the fits for £)„ are shown in Table |2] 

As a first illustration of the level of improvement that can be achieved for the 
energy, Fig. [3] shows the difi^erence between the approximate spectrum with various 
levels of improvement -Bapprox and the exact spectrum inexact, as a function of 77^, 
through the quantity logio(| Ai;!), where l^E = (-Eapprox--Eexact)/^oxact- As expected, 
with each new parameter a new eigenvalue is reproduced, with the concomitant 
reduction in the error. 
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Figure 3. (Color online) Logarithmic plot of the difference between the exact 
spectrum of the two-body problem Eq. ( |88[ l at unitarity, and the approximate 
spectrum obtained with various levels of improvement Nq = 1 — 5, relative to 
the exact spectrum (see text for details), as a function of = Bj/Po' where 
Pq = 2it/L. The original data are discrete; the lines are intended as a guide to the 
eyes. These results correspond to a lattice of side A^^, = 20 and temporal spacing 
6^ = 0.05. 



4-3.2. Contact. The same ideas can be applied to the calculation of Tan's "contact" 
parameter. One of the many ways to define the contact C (see Refs. [5^ [501 [CTl 
is through the derivative of the energy with respect to the inverse scattering length: 

The generalization of the above to finite temperature in the grand canonical 
ensemble is 

on \ 1 /aiogZ\ h"^ 



where is the grand thermodynamic potential, and ^ is the chemical potential. In 
the large- /3 limit, 

C(/3) Co + bc^r' + hc^e-^' (97) 

p— >oo 

where Cq is the ground-state contact, 
47rM51ogZo 

li2 da=^ 



AnMZ.fdE, dE, 



h? \da-^ da- 

Since E may be obtained directly from the logarithm of the partition function, we are 
again in a situation where we require a derivative of the transfer matrix with respect 
to a parameter, in this case . 

In many-body lattice calculations, using these definitions involves the following 
expression: 

^-e'^C/,-.e'^ (100) 
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Table 3. Results of fitting the coeflncients _F„ to the low-energy spectrum of the 
two-body problem at resonance, in a box of side = 16, for an imaginary time 
step = 0.05, in lattice units. 

Nq f, f, f, f. 



1 0.36773 

2 0.14532 0.07568 

3 0.11370 0.02220 

4 0.09659 0.01695 

5 0.08205 0.01278 



0.01957 

0.00415 0.00538 

0.00406 -0.00023 0.00180 



where 



Ua-i = ^— :i sin ct)i. (101) 



The first step towards using these expressions in combination with the 
improvement procedure is to take a formal derivative of T with respect to a^^ in 
the two-particle space, which we can treat exactly: 

In order to compute the change in the exact two-particle energy E2 due to a small 
change in the inverse scattering length, we use the fact that the energies are implicitly 



defined as solutions of Eq. (|43| , which implies 

(103) 



dE^ _ 47r3 dS 



where ry^ = £'2^^/(2"')^: and the derivative on the right-hand side is to be evaluated 
at the corresponding solution of Eq. ( 43 ) . Table [t] in the Appendix shows the first 30 



roots of S{rj) and the corresponding values of dS/drj^. In the derivation of Eq. (103), 
we have assumed that all the effective-range parameters other than the scattering 
length are kept constant. 

Having the exact target spectrum, we proceed by finding the corresponding 



expression in terms of the HS function A{p), which we obtain using Eq. (41) and 
the Feynman-Hellmann theorem: 

^^ = (i?|e -^2V^^e ^Hi^), (104) 
where, as before, we expand in terms of our chosen set of operators, 

^-E^,A.(P.), (105) 

n=0 

and we determine the coefficients F^ by fitting the diagonal matrix elements in the 
right-hand side of Eq. (104 1 to the exact spectrum of Eqs. (102 1 and ( |103[ ). As with 



the energy, the fitting procedure can be reduced to solving a set of linear equations of 
order Nq x Nq . Illustrative results of such a fit are shown in Table [3] 

As with the energy, a first glimpse at the level of improvement that can be 
obtained at this point. This is shown in Fig. |4j where we display the difference AC 
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between the spectrum with various levels of improvement and the exact spectrum, 
divided by the latter. As in Fig. |3j each new parameter allows one to fit a 
new eigenvalue to high accuracy, matching the desired physics beyond the lowest 
momentum. 

Unlike in Fig. |3j the improvement for eigenvalues beyond those explicitly fit is 
limited, breaking down after the eighth or ninth eigenvalue. From that point on toward 
higher energies, little improvement, if any, is observed as the order of the expansion is 
increased. This behavior is only unexpected in the light of Fig. |3j where the situation 
is (surprisingly) much more favorable. 




Figure 4. (Color online) Logarithmic plot of the difference between the exact 
two-body spectrum Eq. ( |102| l at unitarity, and the spectrum obtained with various 
levels of improvement A^q = 1 — 5, relative to the exact spectrum (see text for 
details), as a function of •q'^ = iJj/Po' where Pq = 2n/L. The original data are 
discrete; the lines are intended as a guide to the eyes. These results correspond 
to a lattice of side Nx = 16 and temporal spacing = 0.05. 



Lattice methods for many-body systems 



33 



5. Extracting the ground state 

5.1. Signal-to-noise 

Nonrelativistic systems offer an advantage over relativistic systems concerning tlie 
large Euclidean time limit: because there exist no backward propagating states, the 
entire time extent may be used to recover the ground state. However, there remains a 
serious issue when taking the large Euclidean time limit that goes beyond the simple 
linear algorithmic scaling with /3. This issue will be referred to as the signal-to-noise 
problem. 

In general, one finds that signal-to-noise problems occur for theories which have 
multiparticle states for which the energy per particle is lower than for the states 
one is trying to study. As an example, let us consider a correlation function for a 
single particle in a two-flavor theory with an attractive interaction between different 
flavors, a trivial example for which Monte Carlo calculation is not necessary, but 
which nevertheless illustrates the signal-to-noise problem quite clearly. To determine 
how many conflgurations are necessary to achieve a desired accuracy, it is useful to 
study the signal-to-noise ratio (SNR) as a function of time. To do so, we will follow 
the prescription of Lepage for determining the SNR for nucleon correlators in Lattice 
QCD calculations [161] . The signal is given by 

Ci(/?) = (i^i(0,/3)) 

j Vcl>V[4>]Kii.^.P). (106) 

where 

K^{P)^{l\Bp---B2B^\l) (107) 

and |1) is a single particle state of interest. As shown in Sec. |4j the correlator is 
dominated by the ground state at large Euclidean time, 

(Xi (</.,/?)) 4^)e-K^', (108) 

where z'^\ E^^^ are, respectively, the operator overlap with and energy of the ground- 
state of the n-particle system. For a single particle, we simply have E^''^ = 0, such 
that at large times the correlator will approach a constant. 

Assuming that our sample size is sufficiently large that the central limit theorem 
holds, the noise may be estimated by calculating the standard deviation, 

4(/3)= {{\K,{d,,P)\^)-\{KM.P))?) 

^ ^ V^V[mi{<t>,P)? -\Cm?] ■ (109) 



The first term in Eq. ( 109 1 is the ensemble average of a product of two single- 
particle propagators. This may be recognized as a correlator for two degenerate, 
distinguishable particles. At large Euclidean time, this quantity will be dominated by 
the ground-state of the two-particle system, 

(|ifi(0,/3)P) zl^\-P^o\ (110) 
At zero temperature, two particles in a box with an attractive interacti on wi ll form 



a bound state. This means that Eq < 0, and the first term in Eq. ( 109 ) grows 



exponentially with time. The second term in Eq. ( 109 ) is simply the square of 
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Eq. (106), which is time- independent, 
which dominates. 

We can now form the SNR as 



Thus, in the large-/? hmit it is the first term 



n = 



cm 



N 



cfg 



N. 



cfg 



.(2) 



(111) 



Because Eq < 0, the SNR is exponentially suppressed as a function of Euclidean 
time. Thus, an exponentially large number of configurations will be necessary to see a 
signal at the large times required to extract the ground state from the correlator. An 
example of an effective mass plot displaying a signal-to-noise problem is illustrated in 
Fig. (left panel). 
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Figure 5. Eflfective mass plots displaying two common statistical issues: signal- 
to-noise problem (left), and overlap problem (right). The signal-to-noise problem 
is seen as an exponential growth of error bars with Euclidean time, while an 
overlap problem is seen as a drift with no plateau at large Euclidean time. In 
both cases the ground-state energy is shown as a red line. 

One possibility for redu cing the signal-to-noise problem lies in the time- 
independent prefactor in Eq. (111). The overlap factors Zq"^ depend on how well 



we are able to model the ground-state wavefunction for our source. Thus, the ability 
to create sources which have large overlap with the ground state of the observable 
of interest and poor overlap with the ground state of the system corresponding to 
the variance may greatly increase the number of time steps over which a signal 
can be extracted. This has been documented for the case of nucleon correlators in 
QCD and gives rise to a "golden window" in Euclidean time in which the nucleon 
is in its ground state but exponential degradation of the signal has not yet set in 
|1621ll63Tll64j . In many cases the prefactor is naturally large, for example, in the case 
of two-component unitary fermions the signal corresponds to a scattering state, while 
the noise is associated with deeply bound states [1651 1166] . The natural choice for a 
source will correspond to a scattering state and thus Zq ' , which measures the overlap 
of the scattering state source with the bound states in the variance, will be much less 
than 1. 

Note that this argument assumes that the central limit theorem holds so that 
we may identify the noise as coming from the standard deviation of the operator. 
For many-body operators, this assumption may not in practice be true. Because 
higher-order moments are sensitive to the energies of higher iV-body bound states, 
it is possible for these moments to also grow exponentially with time. For example. 
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the third moment of the probabihty distribution for the single-particle case discussed 
above will be controlled by the energy of the corresponding three-body statc[|] Thus, 
if £;^"' > 3/2E'^ ', the distribut ion will have non-zero skew which grows exponentially 
with time. According to the Berry-Esseen theorem [1671 I168j . the number of 
configurations necessary for the central limit theorem to apply at a given Euclidean 
time becomes, iVcfg > e'^(^-^o^'~3-^o^'). 

The interpretation of such an issue is that the chosen field distribution is not 
peaked for configurations on which the operator is peaked. We will refer to this issue 
as a distribution overlap problem. This problem is particularly worrisome because the 
naive mean that one calculates may differ significantly from the true quantity, though 



the error bars may be deceptively small (see Fig. 5.1 right panel). Such a problem 



will manifest itself in error bars which do not scale with y-^^fg, and in mean values 

which drift outside the error bars as the sample size is changed. 

The traditional approach to solving this problem is to invoke importance 
sampling, in which the statistical overlap of V [(p] with the desired correlation function 
is increased by using the correlation function itself in a reweighting procedure. For 



example, one may work directly with the ratio used for the effective mass, Eq. 81 and 
rewrite it in the following way: 

C^(/3 + A/3) _ /W[0]O(0,/3) ^^^2) 



where 



CnW) !v4>v[(t>] 

CAr((/.,/3 + A/3) 



©(</), /?) = 



Cw(0,/3) 

V[4>] -T'i^lCw (</),/?). (113) 

The r.h.s. of this equation may be directly sampled using a Monte Carlo algorithm 
of one's choice according to the new probability measure V , using the observable C. 
Note that the probability measure is operator-dependent, thus new configurations 
must be generated for each new type of operator one wishes to study. Despite the 
added computational cost and complexity of implementing importance sampling, this 
technique has been used quite successfully in a variety of calculations '90] . 

An alternative to reweighting is to use a better estimator for C(/3) that is free of 
the distribution overlap problem. One such solution is to use the so-called cumulant 
expansion for the energy, which relies on the observation that the distribution of the 
logarithm of the correlator is nearly normally distributed |113j . Such distributions 
have been shown to be ubiquitous in lattice calculations due to the fact that the 
correlator consists of many products of matrices of random numbers [1131 11651 1169) , 
so that for large j3 the distribution becomes very similar in feature to a log-normal 
distribution. Furthermore, the log-normal distribution has been shown to be related to 
systems in which an Efimov effect occurs for distinguishable particles, such as systems 
of unitary fermions [166| . 

§ Due to the lack of (anti-)symmetrization between the fields, the three-body state corresponding 
to the third moment will be composed of particles having three different flavors, with an attractive 
interaction between all particles. In general, the partition function may correspond to a different 
number of flavors from the observable; this is referred to in the lattice QCD literature as a "partially 
quenched" theory. 
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One may expand the average of a generic correlator 



ln(C(/3,0))=^- 



(114) 



where k„ is defined as the n-th cumulant of the distribution for hiC(/3,0). This 
may be recognized as the cumulant generating function evaluated at unity. Note 
that this expansion does not place any assumptions on the particular form of the 
distribution other than that the expansion exists. However, in general an infinite 
number of cumulants need to be included for the equality to hold. 

On the other hand, if In C {jS , (jj) is nearly normally distributed, all cumulants 
higher than will be small and the series may be truncated at a finite order, n^^^. 
Thus the problem has been translated from directly estimating the mean of a quantity 
having a long-tailed distribution, which is exponentially difficult according to the 
Berry-Esseen theorem, to that of estimating small moments of a nearly Gaussian 
distribution. Care must be taken, however, in choosing the order of truncation, and 
any systematic errors that arise due to truncation must be estimated. An example of 
results obtained using the cumulant expansion, including evidence for convergence of 
the series, is shown in Fig. 
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Figure 6. Energy for 50 unitary fermions in a harmonic trap of frequ ency lu, 10 
configurations; fits performed using the cumulant expansion (Eq. |114| l) truncated 
at order n^^^. Inset: conventional effective mass "^^.^{t) = logC(/3)/C(/3 + 1) 
(green) and corresponding fitted cumulant effective mass with n^^^ = 4 (blue). 
Figure from |165| . 



5.2. Interpolating fields 

The requirement for eliminating excited-state contributions may be approximated as 



In 



Zi 



For many-body systems, the spacing between energy levels may become arbitrarily 
small as both collective effects and single quasi-particle excitations require relatively 
small energy to excite. Thus, the number of time steps required to reach the ground 
state may become quite large for 1. Adding to the complication is the fact 
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that larger systems tend to probe higher energies than small systems, so that a very 
fine temporal lattice spacing is required to reduce discretization errors. Finally, the 
exponentially growing signal-to-noise problem previously discussed may severely limit 
the time extents for which a signal may be obtained. 

To remedy this situation, one may construct improved interpolating fields to 
try to maximize the overlap with the ground state. For large systems, the ground- 
state wavefunction becomes increasingly complicated and difficult to parametrize, 
potentially leading to exponentially poor overlap with a chosen interpolating field 
as a function of the number of particles. To illustrate this point, we consider the 
simplest choice for an interpolating field, which is simply a Slater determinant of non- 
interacting single-particle states. Suppose the ground state of the system consists of 
non-interacting two-body states. If the overlap of the two-body ground state with the 
corresponding two-body Slater determinant is e, then the overlap of the iV-body state 
with the iV-body Slater determinant will be , giving exponentially poor overlap 
with the true ground state. 

Two-particle correlations may be built in using an A^/2-body Slater determinant 
of two-body ground states. For the unitary case, the two-body wavefunction is known 
and we may build two-particle correlations using this wavefunction. One possibility is 
to use correlators of the form 

Cjv,.iv,(r) = (det5^^(T)) , (116) 

where 

Sf^T) = (*|i^-i(r,0)®if-i(T,0)|at«J) , (117) 

l^*) is the two-body ground state, and |a"^a^) — ja"^) (E) \a^) are the non-interacting 
single-particle states. This form of the Slater determinant, in which the pairing 
correlations are built into the source and not the sink, can be chosen to avoid the 
order-y increase in computation time required to fully anti-symmetrize both the 
source and sink. It has been shown f l70|, [87] that such a form for the correlation 
functions greatly improves the overlap for systems of unitary fermions. For large, 
fully interacting systems, 3- and higher-body correlations will also contribute to the 
ground-state. These may also be built into the interpolating field, however, each of 
these will increase the computational scaling with V. 

Another way to increase the overlap of the interpolating field onto the ground- 
state wavefunction for systems which possess rotational invariance in the continuum 
is to project the field onto the appropriate irreducible representation of the cubic 
group. By doing so, one may eliminate contributions from low-lying excited states 
with different angular momenta. However, because rotational invariance is broken by 
the cubic lattice, different angular momentum states can mix. For example, if the 
ground state of interest has positive parity and zero total angular momentum, one 
may project onto the representation, which will give overlap with not only the 
J = state, but also J — 4,6,.... Appropriate combinations of sources projected 
onto different representations may be formed to further eliminate contributions from 
undesired angular momenta. For details on creating projection operators see App. [9] 
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6. Systematic effects 



In addition to the statistical errors discussed in the previous section, systematic errors 
arise from the discretization of space and time, as well as from the finite volume of 
space considered. Furthermore, there may be excited state contributions arising from 
the finite temporal extent. In this section, we will largely focus on errors due to spatial 
effects. Temporal discretization errors are subleading to spatial discretization errors, 
br — b'^/M, and may be decreased by increasing the anisotropy parameter M. In 
addition, as discussed below we may improve the transfer matrix to systematically 
remove discretization effects. Temporal discretization effects are also reduced by this 
procedure as the transfer matrix becomes increasingly "perfect" (see Fig. [7| . Finally, 
regarding the issue of finite temporal extent, many improvements to the fitting process 
have been introduced in the lattice QCD literature, such as multi-exponential fits and 
removing excited states by calculating correlation functions using multiple sets of 
sources/sinks [m [TTZl EHl [ISll IM [M] ■ In general, so long as Eq. ( [TTsl ) holds, we 
may be confident that these effects are exponentially suppressed. 
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Figure 7. Ground state energy of TV = 2 + 1 unitary fermions in units of 
Efrcc = {2tt/L)^/M in the A'^ irrep as a function of 1/M for one (left) and four 
(right) tuned couplings with L = 8. Dashed green line corresponds to a fourth 
order polynomial fit in 1/M. Figure from 1891 . 



6.1. Unitary fermions 

To begin our discussion of systematic effects due to spatial discretization and finite 
volume, we will consider the system of two-component unitary fermions. The lack 
of physical scales in this system greatly simplifies the analysis. For a review of the 
introduction of new scales via lattice effective field theory see Ref. [77] . 

Recall that the relevant scale for this case is the Fermi momentum, kp = 
(STT^n)^/'^, where n — N/L^f^^^ is the density. Thus, for a given number of particles 
the size of the box, Lp^yg — bgL, should be considered a fixed physical scale, and the 
lattice spacing, or number of lattice points L, will control spatial discretization errors. 
To reach the continuum limit we require L 3> l/Qj^kp). 

Finite volume errors of the form a/L may be eliminated by tuning the scattering 
length according to the Liischer formalism, which relates the energy levels of two 
particles in a finite box to the desired scattering phase shift. Note that the Liischer 
formalism requires L ^ r^/bg. As discussed below, a lattice theory will intrinsically 
induce rg ^ bg, placing some constraint on the volumes which may be considered. 
However, corrections to the Liischer method due to finite range are exponentially 
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suppressed, while the non-zero effective range also induces errors which obey power- 
law behavior to the energies of unitary fermions, so that in general these corrections 
may be ignored. Finally, we will not discuss the issue of finite size effects in the 
approach to the thermodynamic limit as these are not particular to the discussion of 
lattice methods. 

Discretization errors arise from two sources. The first is due to the discretization 
of the kinetic operator. As discussed in Sec. [2j one may improve this discretization 
by attempting to reproduce the correct single-particle dispersion relation, either by 
adding higher-order derivative terms and matching the energies order-by-order, or by 
calculating the kinetic operator directly in momentum space, allowing one to reproduce 
the exact continuum dispersion relation to all orders up to a momentum cutoff. 

The second source of discretization errors comes from the finite range of the 
interaction induced by the finite lattice spacing. These types of errors may be 
reduced by improving the two-particle sector. One method for doing so is to introduce 
momentum-dependent interactions tuned according to Liischer's formula, as discussed 
at the end of Sec. [2] 

By tuning higher-order energy levels in the two-body spectrum one is also 
effectively tuning higher-order terms in the effective range expansion of p cot 6 to 
zero. This can be verified by plugging the resulting lattice energy eigenvalues back 
into Liischer's formula to recover the effective pcot S seen by the two-particle system. 
This is shown in Fig. [T] where for higher numbers of tuned coefficients p cot 5 becomes 
progressively closer to zero for all momenta. 

Another way to verify this is by looking at the L-dependence of the energy levels 
not tuned to the continuum values [HTj . The expected scaling of these levels, assuming 
the leading term in the effective range expansion has been tuned to 

7rp/6,cot5~7rl/6,r„_ip2" ^ i(27r)2«+ife2«- V„_i77" , (118) 

can be shown to be 

- l) cx i^-^" , (119) 

where rik,r]l are the fcth energy eigenvalue on the lattice and in the continuum 
respectively. This scaling is confirmed in Fig. [8] which shows the L-dependence of 
the untuned levels rj^ and 779. 

While reducing the effects of unwanted s-wave scattering contributions, 
the introduction of momentum-dependent couplings also induces interactions 
corresponding to higher partial waves. The size of the effects of interactions with 
angular momenta £ > may be deduced by studying the energies of the two- 
body lattice eigenstates classified by their irreducible representation (irrep) under the 
octahedral group O (see Appendix I) . Figure [9] shows the deviation of these energy 
levels 77/7/* — 1, where again r],r]* are the energy eigenvalues on the lattice and in the 
continuum respectively, as a function of 77 for the entire spectrum for unitary fermions 
on an L = 8, 16 lattice with No = 2, 4. 

Using the generalization of Liischer's formula for s-wave scattering, we may 
also determine the scattering phase shifts for the higher partial waves. For p-wave 
scattering, if one assumes tanS^ <C tan^i, one finds (see, for example, |137j ): 

p'cotS,ip)=(^y ^vSiv), (120) 
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Figure 8. L-dependence of energy eigenvalues for two particles in a box after 
tuning Nq parameters. Agreement with Eq. |119| for the L-dependence of levels 
?75 and rjg (which were not tuned) shows successful tuning of effective range 
parameters. Figure from 1971 . 
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Figure 9. (Color online) Energy spectrum, given by — log A = jj (^)'^r), of two 
unitary fermions of mass M = 5 for two and four tuned couplings in a finite box of 
size L = 8, 16 (top and bottom, respectively); 7]* are the corresponding eigenvalues 
of the noninteracting theory. Eigenvalues are labeled by dimensionality of the irrep 
to which they belong: square (A), triangle (E), and circle (T), as well as color 
coded according to the lowest orbital angular momentum component contained 
in each irrep: orange (£ = 1), black {£ = 2), cyan {£ = 3), magenta {£ = 4), red 
(£ = 5), green (^ = 6), and yellow (£ = 9). Plots from [89]. 
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10 as a function of rj for 



where 77 and >S'(77) are defined in Eq. 43 Plugging the lattice eigen values obtained for 
the irrep shown in Fig. [9] into the right-hand side of Eq 
prediction for cot Si . This phase shift is plotted in Fig 
L = 16, and for No = 1, 2, 3, 4 and 5. 

To determine the size of the errors for the A'^-body system, we may construct 
the Symanzik effective action [1711 11721 [1731 1174[ 15^ . This is done by considering all 
possible operators formed using the fields and temporal or spatial derivative operators, 
multiplied by unknown, dimensionless coefiicients and appropriate powers of the lattice 
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Figure 10. (Color online) 5i as a function of the dimensionless parameter rj for 
unitary fermions of mass A/ = 5 on an L = 16 lattice obtained from Eq. |120[ 
Figure from |89| . 



spacings to restore the proper dimension. The first operator to appear, assuming the 
leading two-body s-wave contributions have been tuned away as discussed above, will 
be the two-body p-wave interaction discussed above, scaling as L~^. That this scaling 



behavior dominates in the three-body system is verified in Fig. 14 In principle, one 



may tune away such p-wave contributions in the same way as the s-wave contributions, 
by introducing p-wave operators and tuning the coefficients according to the Liischer 
formalism. 

The next class of operators to appear corresponds to iV-body interactions, the 
lowest being the zero-derivative three-body operator. Such an operator will not 
obey a naive scaling law due to the conformal nature of the system, but may be 
shown to scale as L~3.55 ^^^j^ L~^'^^ for £ = 1,0, respectively, using an operator-state 
correspondence |1751 152| . Additional operators corresponding to the four-derivative 
p-wave and c?-wave two-body operators come in at order L~^. 

Knowledge of the scaling of systematic errors with L are not only useful in 
estimating the size of systematic errors, but may also be used in an extrapolation to the 
continuum limit, as seen in Figs. |14|17| For iV = 3, an extrapolation in proved 



to be sufficient, while higher order contributions were negligible. For iV = 4 it was 
found that the contribution from the j>-wave three-body operator was non-negligible, 
so an extrapolation was performed based on the function cq + ci/L^ + C2/L^'^^ ■ 



6.2. External potentials 

The introduction of an external potential, such as a harmonic trap, may also add new 
scales to the problem. Furthermore, the way the potential is included in the transfer 
matrix affects the size of temporal discretization errors. As a simple example, one may 
use Trotter's product formula in the absence of interactions to eliminate the leading 
order contribution, 

T,,t = e-^^/2g-^c/g-^if/2 (^21) 

When interactions are turned on, the order at which temporal discretization errors 
appear depends on how well the interaction has been tuned to remove discretization 
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errors. If we consider the transfer matrix for unitary fermions without an external 
potential, 



/ unitary — ^ ^ 

then we may write the full transfer matrix as 



-b^K/2-b^U/2-b^V,„t^^b^U/2~b^K/2 



(122) 



(123) 



If discretization errors in Ji have been completely eliminated by improvement, one 
may show that temporal discretization errors again come in at 0{b^) [57]- The 
remaining error from temporal discretization may be controlled by tuning the energy 
scale introduced by the potential of interest; for the case of a harmonic potential, this 
scale is the trap frequency, ui. 

For the simple harmonic oscillator (SHO), a new length scale is introduced, 
Lq = (Mw)^/"*. Spatial discretization errors are now controlled by the dimensionless 
ratio Lo/bs, while finite volume errors are controlled by the independent combination 
L/Lo. 

In the continuum limit, finite volume errors for the noninteracting system may be 
computed analytically, since the SHO is separable. A plot of the energy dependence 



of the SHO on L/Lq is shown in Fig. 11 for several low energy single fermion states; at 
large L/Lq, the energies in units of w are just an integer plus the zero point energy (3/2 
for a three-dimensional SHO). For very small volumes, the harmonic potential plays no 
role and the system is effectively a free particle in a finite box, with energie s in creasing 

with decreasing L/Lq. The dashed lines in Fig. 
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proportional to 
this limiting behavior for several SHO states. 
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Figure 11. L/Lq dependence of the SHO energies i^n 

(n = '}2j ) corresponding to tlie single fermion states n = 
(0,0,0),{1,0,0),(1,1,0),(1,1, 1),(2, 1,0),(3, 1, 1), and (4,2,0), wliere 
Lq = (Moj)^/'*. Solid lines indicate an exact continuum limit calculation, 
whereas the data-points indicate simulation results for u) = 0.005 and Lq > 2. 
Dotted lines correspond to free fermions in a finite box (small L/Lq limit). 
Figure from I97| . 



When interactions are turned on, in general the errors as a function of L, Lq 
must be explored numerically. Benchmark results, which will be presented in the 
next subsection, are useful in determining values for bs/Lo and L/Lq for which finite 
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volume and discretization errors are minimal. An example of such a study is shown 
in Fig. [121 
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Figure 12. Ground state energies (in units of uj) as a function of Lg/bs at fixed 
L/ba = 48 for various values of TV, where Lq = (Muj)^^*. Dotted lines are results 
from 11761. Figure from 1971. 



6.3. Comparison with benchmark results 

From testing new code to estimating the size of systematic errors for a given parameter 
set, benchmark results, whether calculated using lattice or other methods such 
as GFMC, are indispensable when developing new methods. Due to the reduced 
computational cost and complexity, few-body systems often serve in this role. Such 
results exist for unitary fermions in both the homogeneous and trapped cases. 
Verification of these results as the variety of methods grows serves to both strengthen 
and constrain the results and give validity to new methods. 

The simplest system to investigate is that for three unitary fermions in a box. 
The energies for the first few zero-momentum states were calculated to high accuracy 
on very large lattices using an unimproved four-fermion interaction in |177j . along 



with a linear extrapolation to the continuum limit (Fig. 13). This measurement was 
confirmed in [971 189] using an improved interaction with four tuned coefficients. An 
extrapolation based on the L^^ scaling of the untuned p-wave operators was performed 



to eliminate the leading contribution to the remaining systematic error (Fig. 14). 

Additional benchmarks for three-body observables, the fermion-dimer scattering 
length and effective range, have been calculated using lattice methods in combination 



with the Liischer formalism (see Sec 
to the dimer binding energy [142 



4]), modified to account for topological corrections 
1179) . The results found using two different 
lattice Hamiltonians with a linear extrapolation to the continuum limit are shown 
in Fig. [15] Also shown is the data computed without correcting for the topological 
volume factor, T{ri); it is clear that the topological factor has a sizeable impact, 
particularly on the result for the effective range. Because the scattering amplitude 
can be reduced to a two-body problem, it can be solved nearly exactly through 
an integral equation formulation, known as the Skorniakov-Ter-Martirosian (STM) 
integral equation |180[ 11811 I182j . Results from the STM equation are included in 



Lattice methods for many-body systems 



44 




Figure 13. Ground (bottom) and first excited (top) state energies for A'^ = 2 + 1 
unitary fermions in a box with zero total angular momentum as a function of b^/L 
as calculated in 11771 . Energies are in units of Eg = 47r^/(2ML^). Lines represent 
a linear fit to the data for b^/L < 1/15, and stars show the energies predicted by 
the Bethe-Peierls model [TtII . 
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Figure 14. Energy of Af = 2 + 1 unitary fermions in a zero total momentum 
eigenstate in units of -EfVoc = 4:Tt^ / {M L^) as a function of 1/L^. Blue data points 
and associated error bars were obtained from numerical simulation, short blue 
dotted lines at L = 8 and L = 10 indicate results from exact diagonalization of 
the three fermion transfer matrix. Red error band indicates the infinite volume 
extrapolation result reported in 1971 189 | using simulation data. Black dashed line 
indicates the infinite volume result of Pricoupenko and Castin reported in |177| . 
Figure from | 97 |. 



Fig. [TSj showing excellent agreement with the lattice results. 

The next system to consider is that for four unpolarized unitary fermions in 
a box. The authors of |183j performed an extensive analysis of this system using 
four different techniques: hamiltonian lattice formalism using iterated eigenvector 
methods, Euclidean lattice formalism with auxiliary-field projection Monte Carlo, 
and continuum diffusion Monte Carlo with fixed and released nodes. The results 
of polynomial extrapolations to the continuum for the lattice methods are shown in 
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Figure 15. Lattice results and continuum extrapolation with statistical 
error bands for the fermion- dimer scattering length (left) and effective range 
(right). All units are given by the binding momentum of the dimer, k. The full 
calculation accounts for the topological volume factor, T{r]) ^ 1. For comparison 
the continuum result obtain via the Skorniakov-Ter-Martirosian equation is also 
presented. Figure from |179) . 



Fig. |16[ The agreement between four different methods gives a very robust benchmark. 
Additionally, this result was verified in [S9 using both an unimproved interaction 
in conjunction with a linear extrapolation to the continuum limit and an improved 
interaction with an extrapolation based on the untuned terms in the Symanzik action 
discussed above (Fig. [l7| . 

Finally, there exist exact results for three unitary fermions in a trap |1751 1184) , 
as well as high-precision solutions to the Schrodinger equation for trapped systems of 
unitary fermions for up to = 6 [1851 1176j . These results have been used to tune the 
parameters L,Lo for the lattice theory described in [ST] (Fig. 12 1. The final results 
from the lattice theory have been shown to agree with the benchmark results to within 
the approximately 1% error bars (see Fig. 18 1. 
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Figure 16. (Color online)Ground state energy for A'^ = 2 + 2 unitary fermions in 
a box as a function of b/L as calculated in |183| in units of Eq = iir'^ / {M L'^) . The 
results are extracted using three methods: an exact Hamiltonian lattice method 
with two different Hamiltonians (blue, red), and a Euclidean lattice method using 
auxiliary field Monte Carlo (green). Lines represent extrapolations to the infinite 
volume limit using a third-order polynomial fit (blue, red) and a linear fit (green), 
with the extrapolated results shown as stars. 
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Figure 17. Ground state energy of A'^ = 2 + 2 unitary fermions as a function 
of bg/L in units of -Efroc = 47r^/(ML^). Error bars include statistical and fitting 
systematic errors combined in quadrature. The red and green bands represent 
fit results to Nq = 1 and = 5 data as discussed in the text, with error 
bands reflecting both statistical and systematic errors. Black dotted lines indicate 
the error band obtained from an infinite volume extrapolation of benchmark 
calculations reported in [183j . Figure from [89j . 
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Figure 18. Fit results for the ground state energies (in units of lu) for 
TV < 6 as a function of L/Lq for two volumes, L = 48,64, and three trap sizes, 
Lq = 7.0, 7.5, 8.0. The result of a correlated fit in N, r of all data is shown as a 
solid line, with a red band showing the combined statistical and systematic errors. 
The results from [TTSl [TTSl (N=3) and [TTSl (N=4-6) are given by dashed lines, 
with any associated error bars shown by hatched regions. Figure from 1971 . 
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7. Selected Results 

In this section we present a few selected lattice results for few- and many-fermion 
systems in the unitary limit. This is only a sample of what the methods can accomplish 
and is by no means fully representative of all the efforts in the field. 

7.1. Unitary fermions in a trap 

As a first demonstration of the techniques outlined in this article to many-body 
systems we present results for the ground-state energies of up to = 70 unpolarized 
unitary fermions in a harmonic trap from j97| . Agreement with benchmark results 
for this system using an improved action with four tuned couplings were presented in 
Sec.lel 

In addition to the increased computational time associated with larger systems, 
many of the issues discussed in previous sections become amplified as the number of 
particles is increased. For example, as more particles are added to the trapped system, 
the wavefunction spreads out in coordinate space, increasing sensitivity to the periodic 
boundary at x = L. Therefore, finite volume errors increase and in general one needs 
to perform an infinite volume extrapolation. To aid in this extrapolation, one may 
use the knowledge that asymptotically the wavefunction for unitary fermions in a trap 
has the same behavior as that for free particles in a trap. Thus, finite volume errors 
will fall off as a Gaussian function, and one may perform a fit of the energies using 
the following form: 

f{L/Lo) - i?co (l - Ae-^(^/^«)') , (124) 

where i?oo is the energy at infinite volume, Lq is the characteristic size of the trap, 
and A, B are fitting parameters. An example of such an extrapolation is shown in 
Fig. [19} 
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Figure 19. Volume dependence of the ground state energies (in units of cu) for 
N = 70. The data points indicate the individual results for six values of L/Lq, 
corresponding to two values of Lq = 7.5, 8.0. The infinite volume extrapolation 
is shown as a solid line, while the band represents the combined statistical and 
fitting systematic errors from extrapolations performed for the two values of Lq 
both separately and combined using Eq. (|124[|. Figure from [97j . 
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The final results for the extrapolated energies as a function of N are shown in 
Fig. [20] in units of the corresponding energies for N noninteracting fermions in a trap. 
These results were generated with the same improved action as that for the few-body 
calculations discussed above using the heat-bath algorithm and utilizing the cumulant 
expansion Eq. (114) to control the large overlap problem. The combined statistical, 



fitting systematic, discretization, and finite volume errors are well under the 1% level 
for large N. 

A comparison with two fixed-node calculations, Diffusion Monte Carlo (FN- 
DMC) and Green's Function Monte Carlo (GFMC), are also shown, along with the 
benchmark results for = 4, 6. In addition to having generally lower results than the 
two variational approaches, the lattice data shows marked shell structure, potentially 
indicating that TV = 70 is not sufficient to reach the thermodynamic limit for the 
trapped system. 
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Figure 20. Ground state energies of A'^ trapped unitary fermions in units of the 
corresponding energies of N trapped noninteracting fermions as a function of N. 
For comparison, we show results from GFMC [I86], FN-DMC [US], and CG [176] 
methods. Figure from |97| . 



7.2. Unitary fermions in a box: energy 

Lattice results from ^89] for the ground-state energies of up to = 66 unpolarized 
unitary fermions in a periodic box in units of the corresponding noninteracting energies 
are presented in Fig. [22] As discussed in the Introduction, this ratio of energies is an 
important universal quantity known as the Bcrtsch parameter. There have been many 
calculations of this quantity, both numerical and theoretical, as well as extractions 
from experiment. A plot of the values reported by various groups is shown in Fig. |21[ 
with values and references for the points in Table |4] 

The results presented in Fig.[22]were generated using an improved action with five 
tuned couplings on lattices ranging from L = 10 — 14. The calculation was performed 
using the heat-bath algorithm in combination with the cumulant expansion to control 
the overlap problem. The results for this calculation are shown in Fig. |22| The result 
of a fit to the region A^ > 34 (representing the thermodynamic limit) is shown as the 
last point in Fig. [21] 

Due to increased statistical issues, involving both signal-to-noise and overlap 
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Figure 21. Historical results for tlie Bertsch parameter determined 

experimentally, by analytic calculation, and by numerical simulation. Numerical 
values and citations are tabulated in Table IZ] Figure from |89| . 



problems, with the untrapped system, the error bars for this calculation are larger 
(~ 2%) than those for the trapped system, however, they are still competitive with 
other sophisticated numerical calculations (see Table |4]). In addition, the use of a 
highly improved action shows no variation in the energies with the volume to within 
the quoted error bars. 
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Figure 22. Ground state energy of untrapped unitary fermions in units of 
Epreeip) = 3NEp{p)/5, where Ef{p) = (37r2p)2/3/2M, as a function of N for 
L = 10, 12 and L = 14. Figure from \S9\. 
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Table 4. Historical results for the Bertsch parameter ^ determined experimentally 
(exp.), by numerical simulation (sim.) and by analytic calculation (anal.), along 
with publication (pub.) date. Values obtained variationally are upper bounds, 
and are indicated with an asterisk; simulation results without a quoted error bar 
should be regarded as approximate. 
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7.5. Unitary fermions in a box: contact 

To illustrate the effect of improved operators in realistic lattice calculations, this 
subsection presents the results of ground-state MC calculations of the contact, for 
an unpolarized system at unitarity. 

Results are shown for 80 particles (40 per spin) in a volume of 10^ lattice points, 
for various improvement levels Nq = 1 — 4. In each case, calculations were performed 
for time directions of extent f3ep = 2.0 — 8.0 (corresponding to roughly between 40 
and 200), subsequently extrapolating to the /? -> oo limit. We have taken 6,- = 0.05, 
and a Slater determinant of plane waves as the starting guess for the ground-state 
wavefunction. For each value of (3, we obtained approximately 400 samples of the 
auxiliary field (f), and the statistics was enhanced by a factor of 20 by inserting 
operators at every other time slice. 

Using extrapolation formulae at leading plus next-to leading order, we obtain 
C/{Nkp) = 3.17(3) in the unimproved case, and C/{Nkp) = 3.31(4), 3.34(4) and 
3.30(4) for Nq — 2, 3, 4, respectively (discarding three data points at the lowest values 
of fitp to capture the asymptotic behavior at large /3; the fits are stable against further 
removal of points at low /3e^). While the latter extrapolations clearly overlap when 
taking the uncertainty into account, they do not overlap with the unimproved case, 
which is clearly consistent with what we see in Fig. [23] at large /3. The remaining 
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Figure 23. (Color online) Contact, in units of Nkp, where kp = {3tt'^N/V)^^^ 
is the Fermi momentum, as a function of the extent of the time direction 
f3€p = rN-rtp, for 80 particles and levels of improvement Nq = 1 — 4. Also 
shown are fits to the = 1 and Nq =4 datasets using the first two terms in the 
asymptotic form of the previous section (see text for details). 



systematic errors, likely due to volume effects for the most part, appear to be only 
as large as 3%, if we consider the most recent estimates of the contact in the ground 
state (~ 3.39, see Ref. |220j '). Remarkably such a small volume as 10'^ already yields 
a result that is close to the best current estimate. 
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8. Summary 

In this short review article we have covered some of the most recent progress on the 
technical aspects of calculating the ground state properties of strongly interacting 
Fermi gases. We have focused on non-rclativistic Hamiltonians with two-body 
interactions, as relevant for ultracold atom experiments and dilute neutron matter 
among other systems. 

Once the problem is placed on a space-time lattice, the partition function is 
written as a path integral over the HS auxiliary field. In order to compute expectation 
values of operators in a non-perturbative fashion, one resorts to stochastic methods 
to evaluate path integrals, assuming that a useful positive semi-definite probability 
can be defined. We have discussed various sampling algorithms to accomplish this, 
from the simplest approaches based on local updates of the auxiliary field, to more 
sophisticated HMC methods, and a very recent technique based on heat-bath-like 
ideas. 

To reduce lattice-spacing effects one must take the low-density limit and/or 

implement improved actions, along the lines of similar efforts in Lattice QCD. We 
have explained how this can be accomplished by designing improved transfer matrices. 
Furthermore, we have shown that that approach can be extended to the construction 
of improved operators. 

We have also presented a small set of illustrative results. 

Finally, the Appendix contains an abridged discussion of some elementary aspects 
of point group theory, as relevant for the discussion of angular momentum (or its 
remnants) on the lattice. 
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9. Appendix I: Group theory on the lattice 

9.1. Introduction and basic definitions 

When placing a field theory on the lattice, i.e. when replacing the spacetime continuum 
with a mesh, one obviously loses a fair amount of symmetry. This is the price to pay 
for going from a problem with infinitely many degrees of freedom to one with a finite 
number. While we have argued that lower densities and improved actions take us closer 
to the continuum limit, we have left for this Appendix the more technical discussion 
of the isometry group associated with the hypercubic lattice. 

In order to understand the degree to which the lattice breaks (or preserves) 
the symmetries of the original theory, it is central to discuss the point group which 
remains as the unbroken group of rotations, namely the hypercubic group Oh, and its 
irreducible representations (irreps) . Indeed, a subspace of states that transform under 
a given irrep should remain invariant under time evolution if the original (continuum) 
Hamiltonian was rotationally invariant. This allows one to classify the scattering 
states according to the discrete equivalent of angular momentum. As we shall see, 
there is only a finite number of irreps of Oh, such that the continuum irreps of the 
rotation group 0(3), characterized by an integer I, are decomposed into direct sums 
of the irreps of Oh- 

Following closely the book by Tinkham [ 226j , we shall first review a series 
of important definitions and theorems (without proof) in the theory of group 
representations. We will bypass the most elementary definitions concerning the groups 
themselves. 

By definition, a representation of a group is a group consisting of concrete 
mathematical entities, such as numbers or matrices. In fact, since we are dealing 
with quantum mechanics, the representations we will discuss are either operators that 
act on an infinite-dimensional Hilbert space (when speaking about the continuum 
limit), or actual matrices (when discussing problems on the lattice). We will have in 
mind the latter case, restricted to square matrices in particular. More specifically, a 
representation associates a matrix T{A) to a given element A of the group in question, 
such that, given any other element B of the group, 

T{A)r{B) = T{AB), (125) 

i.e. the matrices in the representation inherit the multiplication rules from the group. 
It is easy to convince oneself that, in particular, the identity element E must be 
mapped onto the identity matrix, i.e. r(£') = E. It is also easy to see that from a 
given representation T one can obtain another one F' by a similarity transformation 
on all the matrices: F' = STS^^. Those two representations are therefore said to be 
equivalent. 

In general, representations will be reducible, by which we mean that it is possible 
to bring them to a block-diagonal form by a unitary transformation (one and the same 
transformation for the whole original representation), or otherwise irreducible if this 
is not possible. 

In this context, the following lemma is often useful: Any representation 
by matrices with non-vanishing determinant is equivalent through a similarity 
transformation to a representation in terms of unitary matrices. 

This lemma is useful in particular when proving the famous Schur's lemma: Any 
matrix which commutes with all the matrices of an irreducible representation must be a 
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multiple of the identity matrix. If such a non-constant matrix exists, the representation 
must be reducible; if it does not exist, then the representation is irreducible. 

9.2. Orthogonality theorems, characters and conjugacy classes 

The above arc very important, but somewhat formal statements. They allow, however, 
to prove an identity concerning the orthogonality property of the representations that 
is extremely useful in practice: 

^rW(i?);,rO)(i?)„^ = ^s^^s^^, (126) 

R * 

where the sum is over all the group elements, i,j denote particular representations, 
li is the dimension of the i-th representation, and h is the order of the group (i.e. 
the number of elements in it). This equation is valid when i,j denote irreducible, 
inequivalent, unitary representations of the group in question. 

As we shall see, this statement is quite powerful. It can be interpreted 
geometrically if we set i=j and regard R as denoting a coordinate in a vector space 
of dimension h. The above is then a statement of orthogonality between as many 
vectors as matrix elements in the i-ih representation, of which there are If. This is 
valid for all representations, so the above indicates that there are If orthogonal 
vectors. But if the dimension of space is h, then we must have ^ ■ /? < h. As it turns 
out, the equality is valid, so we have 

Y.^^ = h, (127) 

i 

which dramatically constrains the possible dimensions of the collection of irreducible, 
inequivalent, unitary representations of the group. 

We may also obtain a statement about the traces x^^H^) = Trr'^*)(i?) of the 
matrices in the representation. These are usually referred to as characters and they 
have a special name because they are invariant under similarity transformations, such 
that they are a specific signature of each representation. Summing over ji = v and 
a. = j3, we see that characters are orthogonal: 

Y,x'^HRrx'''\R) = h5i,. (128) 

R 

Within a group it is generally possible to identify subsets of elements which, when 
taken together, form a partition (i.e. the subsets are such that each element belongs 
to one and only one subset, and the union of all the subsets yields the full group). 
Each subset is defined by the following (equivalence) relation: Two elements A and B 
in the group are said to be conjugate to each other if an element X exists such that 

A^XBX-'^. (129) 

The subsets in the partition are therefore referred to as conjugacy classes, and contain 
elements connected to each other by the above relation. Clearly, the identity element 
i? is in a class by itself, as it cannot be conjugate to any other element. In the context 
of groups of symmetry transformations, which is the application we have in mind here, 
conjugacy classes roughly correspond to different kinds of operations (e.g. rotations 
about a particular axis). 
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With the above definition, it is clear that the characters of two elements that are 

in the same class will be the same. We may therefore consider the characters of the 
classes Ck within a given representation i instead of those of each of the elements R: 

Y^X^'HRTx^'HR) = Y,x^'\Ckrx^'\Ck)N, = hSij, (130) 

R k 

where A''^. is the number of elements in the fc-th class. Once again, a geometric 
interpretation helps: If we regard the characters as vectors with entries in the space 
of conjugacy classes, then the above indicates that there exist as many orthogonal 
vectors in this space as irreducible representations. But clearly this mimbor cannot 
exceed the dimensionality of the space. Therefore, the number of inequivalent, unitary 
irreps Nirreps is less than or equal to the number of classes N^^. As it turns out, they 
are equal: 

Nirreps = Nq. (131) 

Using the above identities one can also show another orthogonality relation 
between characters: 

Ex^^HCfe)*X«(C;) = ^4;. (132) 

i ^ 

9.3. Decomposition of reducible representations 

If you were wondering when all of the above would become useful, this is it. 

The character of a reducible representation is obviously the sum of the characters 
of all the representations it contains, i.e. 

x{R) = Y.a^X^^\R) (133) 

3 

where is the multiplicity of the j-th irrep, i.e. the number of times it appears. 
The above orthogonality relations can therefore help us identify the irrep content of a 
given representation: 

= /i-' = h-'Y.^^^^'^{Ckrx{Ck). (134) 

R k 

9.4- Projection operators for irreducible representations 

the 

basis element >Pk \ the other ones are usually referred to as partners. By definition, 
the action of an element of the group on ' can be expanded as a linear combination 
in the basis of that function and its partners: 

^'K<^i^^ = E^A^r«(^)^.- (135) 

A=l 

Multiplying by r^'\R)l, summing over R, and using the orthogonality theorems, 
we obtain 

ErW(i?)^,,,p«<^0) = (136) 

R '■j 



Let the basis of the j-ih irrep be denoted by where a = 1, Z^ . Given a particular 
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Clearly then, applying to a function the operator 

■^i^ = ^Er^'n^)U^i^ (137) 

R 

yields zero unless the function has a projection on the K-th basis element of the j-th 
representation, in which case the result is transferred onto the A-th direction. For this 
reason these operators are sometimes referred to as transfer operators. The special 
case X = K yields a projection operator, as in that case we have 

^!^M^=^i^^- (138) 

While the above identities are useful when we already know the explicit form of 
the r matrices, one can still obtain practical results by knowing just the characters. 
Indeed, setting A = k in Eq. 137 and summing over k we obtain a new projection 
operator 

^0) ^^^(^ ^^AJ2^U)^jiyp^ (139) 

K R 

which projects onto the j-th irrep. 



9.5. The octahedral group 

By far the most common choice in lattice field theory calculations is to discretize 
spacetime as a hypercubic mesh. Therefore the point group one is concerned with 
is the so-called octahedral group O. This group contains the proper (no parity 
transformations) rotations which take a cube or an octahedron into itself. In this 
case h = 24, i.e. it has 24 elements, which include, apart from the identity (E), eight 
rotations about cube body diagonals (8C3), nine rotations around the X,Y,Z axes 
(3C2 plus 6C4), and six rotations about axes parallel to face diagonals (6C2). Should 
we choose to include the inversion operation as well, which form the group i, we will 
then obtain the group Of^, where = O x i. This is referred to as the full octahedral 
group, and with 48 elements it is the largest of all the point groups. 



Table 5. Tabic of characters of the octahedral group O. 

E 8C3 3C2 6C2 ~6C^ 



1 1 1 1 1 

^2 1 1 1-1-1 

£; 2 -1 2 

Ti 3 -1-1 1 

T2 3 -1 1 -1 



The table of characters of O is shown in Table [5] The classes appear in the top 
row, and the irreps along the first column. The dimensionality of each irrep can be 
read from the second column. 

While we have only 5 irreps for O, in the continuum we have SO (3), and therefore 
an infinite number of irreps, parametrized by an integer i, with dimensions given by 
2£ + 1 and spanned by the spherical harmonics F/" . It is therefore natural to ask: 
within each irrep of 5*0(3), how is the dimensionality split among the finite number 
of irreps of the octahedral group when we put the problem on the lattice? 
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In order to answer this question, we may simply use the character orthogonahty 
relations on the ^-th representation of 50(3). We need only know the character of the 
elements of O, now interpreted as elements of S0{3), in that representation, which 
are given by 

, s'm(£ + h)a , ^ 

^(") = ^ii^- 

Thus, 



^ = 0,3,... 
f =1,4,... 
f = 2,5,... 

1 ^ = 0,1,4,5,... 
-1 £ = 2,3,6,7,... 

Using these results one may then decompose the s-wave, p-wave, etc., 
representations of 50(3) in terms of the irreps of the octahedral group. The (partial) 
result of such a decomposition is shown in Table |6] A considerably more thorough 
treatment of this problem appears in Ref. |137j . From the table we see that only 
the A-^ representation contributes to s-wave scattering, but it also appears in higher 
waves. For £ > A the various irreps of O will in general appear more than once. 




= x(7r/2) = 



Table 6. Characters and decomposition of irreps of 50(3) into irreps of the 
octahedral group O. 
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Table 7. First 30 roots of 5 (77), and dS/dri^ evaluated at those roots. 
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